Merge pull request #55 from BattModels/concentration_corrections

add ability to model activities as a function of composition and correction for equilibrium voltage

Project.toml

@@ -1,7 +1,7 @@
 name = "ElectrochemicalKinetics"
 uuid = "a2c6e634-85ca-418a-9c67-9b5417ce2d04"
 authors = ["Rachel Kurchin <rkurchin@cmu.edu>", "Holden Parks <hparks@andrew.cmu.edu>", "Dhairya Gandhi <dhairya@juliacomputing.com>"]
-version = "0.1.5"
+version = "0.2.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/fitting.jl

@@ -21,8 +21,16 @@ function janky_log_loss(y, y_pred)
 end
 
 # helper fcn to make sure starting guess for `overpotential` has correct sign
-function _get_guess(k, model)
-    guess = fill(1f-1, max(length(k), length(model)))
+function _get_guess(init_guess, k, model, a_r, a_o)
+    guess = init_guess
+    max_l = max(length(k), length(model), length(a_r), length(a_o))
+    if init_guess isa Real
+        guess = fill(init_guess, max_l)
+    else
+        if length(init_guess) != max_l
+            @warn "Length mismatch between guess and required solution length"
+        end
+    end
     # make sure guess has correct sign
     if k isa Real
         if k == 0
@@ -45,9 +53,10 @@ Given values for current/rate constant and specified model parameters, find the
 
 NOTE that this currently only solves for net reaction rates.
 """
-function overpotential(k, model::KineticModel, guess = _get_guess(k, model); T = 298, loss = janky_log_loss, autodiff = true, verbose=false, warn=true, kwargs...)
+function overpotential(k, model::KineticModel; a_r=1.0, a_o=1.0, guess = 0.1, T = 298, loss = janky_log_loss, autodiff = true, verbose=false, warn=true, kwargs...)
     # wherever k=0 we can shortcut since the answer has to be 0
     k_solve = k
+    guess = _get_guess(guess, k, model, a_r, a_o)
     if k==0 # scalar k=0, possibly vector model
         if verbose
             println("shortcutting full solve")
@@ -58,28 +67,29 @@ function overpotential(k, model::KineticModel, guess = _get_guess(k, model); T =
             println("shortcutting part of solve")
         end
         k_solve = k[k.!=0]
+        guess_solve = guess[k.!=0]
     end
 
     function compare_k!(storage, V)
-        storage .= loss(k, rate_constant(V, model; T = T, kwargs...))
+        storage .= loss(k, rate_constant(V, model; a_r=a_r, a_o=a_o, T = T, kwargs...))
     end
     Vs = if autodiff
         function myfun!(F, J, V)
             # println("V=",V)
             # println("loss=", loss(k, rate_constant(V, model; T=T, kwargs...)))
             if isnothing(J)
-                F .= loss(k_solve, rate_constant(V, model; T = T, kwargs...))
+                F .= loss(k_solve, rate_constant(V, model; a_r=a_r, a_o=a_o, T = T, kwargs...))
             elseif isnothing(F) && !isnothing(J)
                 gs = Zygote.gradient(V) do V
                     Zygote.forwarddiff(V) do V
-                        loss(k_solve, rate_constant(V, model; T = T, kwargs...)) |> sum
+                        loss(k_solve, rate_constant(V, model; a_r=a_r, a_o=a_o, T = T, kwargs...)) |> sum
                     end
                 end[1]
                 J .= diagm(gs)
             else
                 y, back = Zygote.pullback(V) do V
                     Zygote.forwarddiff(V) do V
-                        loss(k_solve, rate_constant(V, model; T = T, kwargs...)) |> sum
+                        loss(k_solve, rate_constant(V, model; a_r=a_r, a_o=a_o, T = T, kwargs...)) |> sum
                     end
                 end
                 F .= y
@@ -101,11 +111,12 @@ function overpotential(k, model::KineticModel, guess = _get_guess(k, model); T =
     else
         sol_return = sol
     end
+    if length(sol_return) == 1
+        sol_return = first(sol_return)
+    end
     sol_return
 end
 
-overpotential(k::Real, model::KineticModel{T}, guess=0.1; kwargs...) where T<:Real = overpotential([k], model, [0.1]; kwargs...)[1]
-
 inv!(x) = x .= inv.(x)
 
 # TODO: test this
@@ -126,7 +137,7 @@ end
 
 
 # multiple models, one k value (used in thermo example)
-overpotential(k::Real, models::Vector{<:KineticModel}, guess=fill(0.1, length(k)); kwargs...) = overpotential.(Ref(k), models, Ref(guess); kwargs...)
+# overpotential(k::Real, models::Vector{<:KineticModel}, guess=fill(0.1, length(k)); kwargs...) = overpotential.(Ref(k), models, Ref(guess); kwargs...)
 
 fitting_params(t::Type{<:NonIntegralModel}) = fieldnames(t)
 fitting_params(::Type{MarcusHushChidsey}) = (:A, :λ)

src/kinetic_models/AsymptoticMarcusHushChidsey.jl

@@ -27,12 +27,5 @@ function rate_constant(V_app, amhc::AsymptoticMarcusHushChidsey, ox::Bool; T = 2
     return kB * T .* amhc.A .* pref .* erfc.(arg)
 end
 
-# direct dispatch for net rates
-function rate_constant(V_app, amhc::AsymptoticMarcusHushChidsey; T = 298)
-    η = V_app / (kB * T)
-    λ_nondim = amhc.λ / (kB * T)
-    a = 1 .+ sqrt.(λ_nondim)
-    arg = (λ_nondim .- sqrt.(a .+ η.^2)) ./ (2 .* sqrt.(λ_nondim))
-    pref = sqrt.(π .* λ_nondim) .* tanh.(η/2)
-    return kB * T * amhc.A .* pref .* erfc.(arg)
-end
+rate_constant(V_app, amhc::AsymptoticMarcusHushChidsey, ::Val{true}; T=298) = rate_constant(V_app, amhc, true; T=T)
+rate_constant(V_app, amhc::AsymptoticMarcusHushChidsey, ::Val{false}; T=298) = rate_constant(V_app, amhc, false; T=T)

src/kinetic_models/kinetic_models.jl

@@ -64,18 +64,21 @@ abstract type IntegralModel{T} <: KineticModel{T} end
 # integrand(km::IntegralModel, V_dl, ox::Bool; kwargs...) =
 #     error("An integral-based kinetic model must dispatch the `integrand` function!")
 
-# TODO: check that this passes through both kT and V_q appropriately
+# TODO: check that this passes through V_q appropriately
 # dispatch for net rates
-integrand(km::IntegralModel, V; kwargs...) =
-    E -> integrand(km, V, Val(true); kwargs...)(E) - integrand(km, V, Val(false); kwargs...)(E)
+function integrand(km::IntegralModel, V; a_r=1.0, a_o=1.0, T=298, kwargs...)
+    eq_V = kB * T .* log.(a_r ./ a_o)
+    E -> a_r .* integrand(km, V .- eq_V, Val(true); kwargs...)(E) .- a_o .* integrand(km, V .- eq_V, Val(false); kwargs...)(E)
+end
+
 integrand(km::IntegralModel, V, ox::Bool; kwargs...) = integrand(km, V, Val(ox); kwargs...)
 
 """
     rate_constant(V_app, model::KineticModel, ox::Bool; kwargs...)
-    rate_constant(V_app, model::KineticModel; kwargs...)
+    rate_constant(V_app, model::KineticModel; a_r=1.0, a_o=1.0, kwargs...)
     rate_constant(E_min, E_max, V_app, model::MarcusHushChidseyDOS, calc_cq::Bool=false; C_dl = 10.0, Vq_min = -0.5, Vq_max = 0.5, kwargs...)
 
-Compute the rate constant k predicted by a given kinetic model at a applied voltage `V_app`. If a flag for reaction direction `ox` is supplied, `true` gives the oxidative and `false` the reductive direction, while omitting this flag will yield net reaction rate (absolute value thereof).
+Compute the rate constant k predicted by a given kinetic model at a applied voltage `V_app`. If a flag for reaction direction `ox` is supplied, `true` gives the oxidative and `false` the reductive direction, while omitting this flag will yield net reaction rate. In the net rate case, activities of reduced and oxidized species (`a_r` and `a_o`, respectively) may also be supplied (default values are unity).
 
 If the model is an `IntegralModel`, integration bounds `E_min` and `E_max` may be supplied as kwargs. Integration is done via GK quadrature.
 
@@ -84,8 +87,10 @@ If calc_cq flag is passed, optionally compute voltage shifts due to quantum capa
 rate_constant(V_app, model::NonIntegralModel, ox::Bool; kwargs...) = rate_constant(V_app, model, Val(ox); kwargs...)
 
 # dispatch for net rates
-rate_constant(V_app, model::NonIntegralModel; kwargs...) =
-    rate_constant(V_app, model, Val(true); kwargs...) - rate_constant(V_app, model, Val(false); kwargs...)
+function rate_constant(V_app, model::NonIntegralModel; a_r=1.0, a_o=1.0, T=298, kwargs...)
+    eq_V = kB * T .* log.(a_r./a_o)
+    a_r .* rate_constant(V_app .- eq_V, model, Val(true); T=T, kwargs...) .- a_o .* rate_constant(V_app .- eq_V, model, Val(false); T=T, kwargs...)
+end
 
 # TODO: add tests that both args and kwargs are correctly captured here (also for the Val thing)
 # "callable" syntax
@@ -97,10 +102,11 @@ function rate_constant(
     args...; # would just be the ox flag, if present
     T = 298,
     E_min = -100 * kB * T,
-    E_max = 100 * kB * T
+    E_max = 100 * kB * T,
+    kwargs... # e.g. a_o, a_r
 )
     n, w = scale(E_min, E_max)
-    f = integrand(model, V_app, args...; T = T)
+    f = integrand(model, V_app, args...; T = T, kwargs...)
     sum(w .* f.(n))
 end
 

src/phase_diagrams.jl

@@ -17,24 +17,35 @@ g_thermo(x; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) = @.
 
 prefactor(x, intercalate::Bool) = intercalate ? (1 .- x) : x
 
+# TODO: add support for other thermodynamic models than ideal mixing...
 """
 µ and g with kinetic constributions, can be modeled using any <:KineticModel object
 
 These functions return single-argument functions (to easily use common-tangent function below while
 still being able to swap out model parameters by calling "function-builders" with different arguments).
+
+Notably, expressions for computing activity as a function of composition `x` for the oxidized and reduced states can be supplied via the activity_function_o and activity_function_r.
 """
-function µ_kinetic(I, km::KineticModel; intercalate=true, warn=true, T=room_T, kwargs...)
+function µ_kinetic(I, km::KineticModel; 
+    activity_function_o=x-> 1 .- x, # default a la Bazant
+    activity_function_r=x-> 1 .- x, 
+    warn=true, 
+    T=room_T, 
+    kwargs...)
     thermo_term(x) = μ_thermo(x; T=T, kwargs...)
-    μ(x::Real) = thermo_term(x) .+ overpotential(I/prefactor(x, intercalate), km, T=T, warn=warn)
-    μ(x::AbstractVector) = thermo_term(x) .+ overpotential(I./prefactor(x, intercalate), km, T=T, warn=warn)
+    μ(x::Real) = thermo_term(x) .+ overpotential(I, km; a_r=activity_function_r(x), a_o=activity_function_o(x), T=T, warn=warn)[1]
+    μ(x::AbstractVector) = thermo_term(x) .+ overpotential(I, km; a_r=activity_function_r(x), a_o=activity_function_o(x), T=T, warn=warn)
     return μ
 end
 
-function g_kinetic(I, km::KineticModel; intercalate=true, warn=true, T=room_T, kwargs...)
+function g_kinetic(I, km::KineticModel; 
+    activity_function_o=x-> 1 .- x, 
+    activity_function_r=x-> 1 .- x,
+    warn=true, T=room_T, kwargs...)
     thermo_term(x) = g_thermo(x; T=T, kwargs...)
     #TODO: gradient of this term is just value of overpotential(x)
     function kinetic_term(x)
-        f(x) = ElectrochemicalKinetics.overpotential.(I ./prefactor(x, intercalate), Ref(km), T=T, warn=warn)
+        f(x) = ElectrochemicalKinetics.overpotential(I, km; a_r=activity_function_r(x), a_o=activity_function_o(x), T=T, warn=warn)
         n, w = ElectrochemicalKinetics.scale_coarse(zero.(x), x)
         map((w, n) -> sum(w .* f(n)), eachcol(w), eachcol(n))
     end
@@ -46,17 +57,17 @@ end
 # TODO: figure out T passthrough issue
 # zeros of this function correspond to pairs of x's satisfying the common tangent condition for a given µ function
 # case where we just feed in two points (x should be of length two)
-function common_tangent(x::Vector, I, km::KineticModel; intercalate=true, kwargs...)
-    g = g_kinetic(I, km; intercalate=intercalate, kwargs...)
-    µ = µ_kinetic(I, km; intercalate=intercalate, kwargs...)
+function common_tangent(x::Vector, I, km::KineticModel; kwargs...)
+    g = g_kinetic(I, km; kwargs...)
+    µ = µ_kinetic(I, km; kwargs...)
     [(g(x[2]) - g(x[1]))/(x[2] - x[1]) .- μ(x[1]), μ(x[2])-μ(x[1])]
 end
 
 # TODO: see if we can speed this up with gradients? And/or if it's even needed for integral cases
-function find_phase_boundaries(I, km::KineticModel; guess=[0.05, 0.95], intercalate=true, verbose=false, kwargs...)
+function find_phase_boundaries(I, km::KineticModel; guess=[0.05, 0.95], verbose=false, kwargs...)
 
     function myct!(storage, x)
-        res = common_tangent(x, I, km; intercalate=intercalate, kwargs...)
+        res = common_tangent(x, I, km; kwargs...)
         storage[1] = res[1]
         storage[2] = res[2]
     end
@@ -73,20 +84,21 @@ NOTE 1: appropriate values of `I_step` depend strongly on the prefactor of your
 
 NOTE 2: at lower temperatures (<=320K or so), ButlerVolmer models with the default thermodynamic parameters have a two-phase region at every current, so setting a finite value of I_max is necessary for this function to finish running.
 """
-function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, intercalate=true, start_guess=[0.05, 0.95], tol=5e-3, kwargs...)
+function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, start_guess=[0.05, 0.95], tol=5e-3, kwargs...)
     I = I_start
-    pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=start_guess, kwargs...)
+    pbs_here = find_phase_boundaries(I, km; guess=start_guess, kwargs...)
     pbs = pbs_here'
     I_vals = [I_start]
     pb_diff = [0.0, 0.0]
-    while abs(pbs_here[2] - pbs_here[1]) > tol && I < I_max
+    # while within tolerance and "on same side" of max
+    while abs(pbs_here[2] - pbs_here[1]) > tol && (I_max - I) * I_step > 0
         I = I + I_step
         if verbose
             println("Solving at I=", I, "...")
         end
         try
             pbs_old = pbs_here
-            pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_old .+ pb_diff, kwargs...)
+            pbs_here = find_phase_boundaries(I, km; guess=pbs_old .+ pb_diff, kwargs...)
             pb_diff = pbs_here .- pbs_old
             # TODO: check that they haven't crossed over
             pbs = vcat(pbs, pbs_here')

test/phase_diagram_tests.jl

@@ -5,7 +5,7 @@ muoB = 0.03
 T = 298
 
 # TODO: add tests where we change Ω, etc.
-# TODO: add tests for deintercalation case
+# TODO: add tests for other activity functions
 
 @testset "Basic free energy functions" begin
     # enthalpy of mixing
@@ -144,7 +144,7 @@ end
             # they should get "narrower" with temperature
             @test v2[1] > v1[1]
             @test v2[2] < v1[2]
-            v3 = find_phase_boundaries(400, km, guess=[0.1,0.8])
+            v3 = find_phase_boundaries(400, km, guess=[0.1,0.7])
             # ...and also with current
             @test v3[1] > v1[1]
             @test v3[2] < v1[2]
@@ -167,10 +167,11 @@ end
         v1 = find_phase_boundaries(100, mhc, guess=v1_vals[:AsymptoticMarcusHushChidsey])
         @test all(isapprox.(v1, [0.04967204036, 0.81822937543], atol=1e-5))
 
+        # values for these last two taken to be correct as of 2022-10-18
         v2 = find_phase_boundaries(100, mhc, T=350, guess=v2_vals[:AsymptoticMarcusHushChidsey])
-        @test all(isapprox.(v2, [0.075615, 0.8171636], atol=1e-5))
+        @test all(isapprox.(v2, [0.09011, 0.77146], atol=1e-5))
 
         v3 = find_phase_boundaries(400, mhc, guess=v3_vals[:AsymptoticMarcusHushChidsey])
-        @test all(isapprox.(v3, [0.1467487, 0.5704125], atol=1e-5))
+        @test all(isapprox.(v3, [0.1467, 0.57035], atol=1e-5))
     end
 end