A bunch of tweaks, plus some fixes to pass tests

- add activity_function_r, activity_function_o parameters to model activities as a function of composition (they default to the Bazant expressions)
- finish renaming c_o, c_r -> a_o, a_r (concentration -> activity)
- remove separate dispatch of `overpotential` for scalar k/model since activities can still be vectors but those are kwargs
- remove net rate dispatch for aMHC since it only works if activities are equal to each other

Signed-off-by: Rachel Kurchin <rkurchin@cmu.edu>

src/fitting.jl

@@ -21,9 +21,9 @@ function janky_log_loss(y, y_pred)
 end
 
 # helper fcn to make sure starting guess for `overpotential` has correct sign
-function _get_guess(init_guess, k, model, c_r, c_o)
+function _get_guess(init_guess, k, model, a_r, a_o)
     guess = init_guess
-    max_l = max(length(k), length(model), length(c_r), length(c_o))
+    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
@@ -53,10 +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; c_r=1.0, c_o=1.0, init_guess = 0.1, 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(init_guess, k, model, c_r, c_o)
+    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")
@@ -71,25 +71,25 @@ function overpotential(k, model::KineticModel; c_r=1.0, c_o=1.0, init_guess = 0.
     end
 
     function compare_k!(storage, V)
-        storage .= loss(k, rate_constant(V, model; c_r=c_r, c_o=c_o, 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; c_r=c_r, c_o=c_o, 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; c_r=c_r, c_o=c_o, 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; c_r=c_r, c_o=c_o, 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
@@ -111,11 +111,12 @@ function overpotential(k, model::KineticModel; c_r=1.0, c_o=1.0, init_guess = 0.
     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; guess=[0.1], kwargs...)[1]
-
 inv!(x) = x .= inv.(x)
 
 # TODO: test this

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

@@ -67,7 +67,7 @@ abstract type IntegralModel{T} <: KineticModel{T} end
 # TODO: check that this passes through V_q appropriately
 # dispatch for net rates
 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)
+    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
 
@@ -89,7 +89,7 @@ rate_constant(V_app, model::NonIntegralModel, ox::Bool; kwargs...) = rate_consta
 # dispatch for net rates
 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); kwargs...) .- a_o .* rate_constant(V_app .- eq_V, model, Val(false); kwargs...)
+    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)
@@ -102,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; 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, km; c_r=x/((1-x)), c_o=1, T=T, warn=warn)[1]
-    μ(x::AbstractVector) = thermo_term(x) .+ overpotential(I, km; c_r=x ./ ((1 .-x)), c_o=1, 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; 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, km; c_r=x ./ ((1 .- x)), c_o=1, 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
@@ -73,7 +84,7 @@ 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; guess=start_guess, kwargs...)
     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