diff --git a/LICENSE b/LICENSE index 3562176..098f087 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright 2021, Rachel Kurchin and Holden Parks +Copyright 2022, Rachel Kurchin, Holden Parks, and Dhairya Gandhi Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: diff --git a/Manifest.toml b/Manifest.toml deleted file mode 100644 index 79930ae..0000000 --- a/Manifest.toml +++ /dev/null @@ -1,1054 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -manifest_format = "2.0" - -[[deps.AbstractFFTs]] -deps = ["ChainRulesCore", "LinearAlgebra"] -git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4" -uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" -version = "1.1.0" - -[[deps.Adapt]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f" -uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.3.3" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" - -[[deps.ArrayInterface]] -deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] -git-tree-sha1 = "d49f55ff9c7ee06930b0f65b1df2bfa811418475" -uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "4.0.4" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[deps.AxisAlgorithms]] -deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] -git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7" -uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" -version = "1.0.1" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[deps.Bzip2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" -uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.8+0" - -[[deps.Cairo_jll]] -deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2" -uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" -version = "1.16.1+1" - -[[deps.ChainRules]] -deps = ["ChainRulesCore", "Compat", "IrrationalConstants", "LinearAlgebra", "Random", "RealDot", "SparseArrays", "Statistics"] -git-tree-sha1 = "098b5eeb1170f569a45f363066b0e405868fc210" -uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2" -version = "1.27.0" - -[[deps.ChainRulesCore]] -deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "c9a6160317d1abe9c44b3beb367fd448117679ca" -uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.13.0" - -[[deps.ChangesOfVariables]] -deps = ["ChainRulesCore", "LinearAlgebra", "Test"] -git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1" -uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" -version = "0.1.2" - -[[deps.ColorSchemes]] -deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"] -git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640" -uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.17.1" - -[[deps.ColorTypes]] -deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" -uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.11.0" - -[[deps.Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] -git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" -uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.8" - -[[deps.CommonSubexpressions]] -deps = ["MacroTools", "Test"] -git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" -uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" -version = "0.3.0" - -[[deps.Compat]] -deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.41.0" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" - -[[deps.Contour]] -deps = ["StaticArrays"] -git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7" -uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" -version = "0.5.7" - -[[deps.DataAPI]] -git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8" -uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.9.0" - -[[deps.DataStructures]] -deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.11" - -[[deps.DataValueInterfaces]] -git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" -uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" -version = "1.0.0" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[deps.DelimitedFiles]] -deps = ["Mmap"] -uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" - -[[deps.DiffResults]] -deps = ["StaticArrays"] -git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805" -uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" -version = "1.0.3" - -[[deps.DiffRules]] -deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0" -uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.10.0" - -[[deps.Distances]] -deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04" -uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" -version = "0.10.7" - -[[deps.Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" - -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.6" - -[[deps.Downloads]] -deps = ["ArgTools", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" - -[[deps.EarCut_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d" -uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5" -version = "2.2.3+0" - -[[deps.Expat_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "ae13fcbc7ab8f16b0856729b050ef0c446aa3492" -uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.4.4+0" - -[[deps.FFMPEG]] -deps = ["FFMPEG_jll"] -git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8" -uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" -version = "0.4.1" - -[[deps.FFMPEG_jll]] -deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] -git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f" -uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" -version = "4.4.0+0" - -[[deps.FastGaussQuadrature]] -deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "58d83dd5a78a36205bdfddb82b1bb67682e64487" -uuid = "442a2c76-b920-505d-bb47-c5924d526838" -version = "0.4.9" - -[[deps.FillArrays]] -deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] -git-tree-sha1 = "4c7d3757f3ecbcb9055870351078552b7d1dbd2d" -uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.13.0" - -[[deps.FiniteDiff]] -deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "ec299fdc8f49ae450807b0cb1d161c6b76fd2b60" -uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.10.1" - -[[deps.FixedPointNumbers]] -deps = ["Statistics"] -git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" -uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" -version = "0.8.4" - -[[deps.Fontconfig_jll]] -deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03" -uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" -version = "2.13.93+0" - -[[deps.Formatting]] -deps = ["Printf"] -git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" -uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" -version = "0.4.2" - -[[deps.ForwardDiff]] -deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "1bd6fc0c344fc0cbee1f42f8d2e7ec8253dda2d2" -uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.25" - -[[deps.FreeType2_jll]] -deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9" -uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" -version = "2.10.4+0" - -[[deps.FriBidi_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91" -uuid = "559328eb-81f9-559d-9380-de523a88c83c" -version = "1.0.10+0" - -[[deps.GLFW_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"] -git-tree-sha1 = "51d2dfe8e590fbd74e7a842cf6d13d8a2f45dc01" -uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" -version = "3.3.6+0" - -[[deps.GR]] -deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "RelocatableFolders", "Serialization", "Sockets", "Test", "UUIDs"] -git-tree-sha1 = "9f836fb62492f4b0f0d3b06f55983f2704ed0883" -uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.64.0" - -[[deps.GR_jll]] -deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "a6c850d77ad5118ad3be4bd188919ce97fffac47" -uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.64.0+0" - -[[deps.GeometryBasics]] -deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] -git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c" -uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" -version = "0.4.1" - -[[deps.Gettext_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" -uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" -version = "0.21.0+0" - -[[deps.Glib_jll]] -deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06" -uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" -version = "2.68.3+2" - -[[deps.Graphite2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" -uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" -version = "1.3.14+0" - -[[deps.Grisu]] -git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" -uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" -version = "1.0.2" - -[[deps.HTTP]] -deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"] -git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a" -uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.9.17" - -[[deps.HarfBuzz_jll]] -deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] -git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3" -uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" -version = "2.8.1+1" - -[[deps.IRTools]] -deps = ["InteractiveUtils", "MacroTools", "Test"] -git-tree-sha1 = "7f43342f8d5fd30ead0ba1b49ab1a3af3b787d24" -uuid = "7869d1d1-7146-5819-86e3-90919afe41df" -version = "0.4.5" - -[[deps.IfElse]] -git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" -uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -version = "0.1.1" - -[[deps.IniFile]] -git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625" -uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" -version = "0.5.1" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.Interpolations]] -deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] -git-tree-sha1 = "b15fc0a95c564ca2e0a7ae12c1f095ca848ceb31" -uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.13.5" - -[[deps.InverseFunctions]] -deps = ["Test"] -git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65" -uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.2" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.1.1" - -[[deps.IterTools]] -git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5" -uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" -version = "1.4.0" - -[[deps.IteratorInterfaceExtensions]] -git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" -uuid = "82899510-4779-5014-852e-03e436cf321d" -version = "1.0.0" - -[[deps.JLLWrappers]] -deps = ["Preferences"] -git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.4.1" - -[[deps.JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" -uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.3" - -[[deps.JpegTurbo_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba" -uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" -version = "2.1.2+0" - -[[deps.LAME_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" -uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" -version = "3.100.1+0" - -[[deps.LZO_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" -uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" -version = "2.10.1+0" - -[[deps.LaTeXStrings]] -git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" -uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" -version = "1.3.0" - -[[deps.Latexify]] -deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"] -git-tree-sha1 = "a6552bfeab40de157a297d84e03ade4b8177677f" -uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.15.12" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" - -[[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.Libffi_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290" -uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" -version = "3.2.2+1" - -[[deps.Libgcrypt_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"] -git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae" -uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" -version = "1.8.7+0" - -[[deps.Libglvnd_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"] -git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf" -uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29" -version = "1.3.0+3" - -[[deps.Libgpg_error_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9" -uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" -version = "1.42.0+0" - -[[deps.Libiconv_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778" -uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.1+1" - -[[deps.Libmount_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73" -uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.35.0+0" - -[[deps.Libtiff_jll]] -deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"] -git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9" -uuid = "89763e89-9b03-5906-acba-b20f662cd828" -version = "4.3.0+0" - -[[deps.Libuuid_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066" -uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.36.0+0" - -[[deps.LineSearches]] -deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] -git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd" -uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -version = "7.1.1" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.LogExpFunctions]] -deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.6" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[deps.MacroTools]] -deps = ["Markdown", "Random"] -git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf" -uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.9" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MbedTLS]] -deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"] -git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe" -uuid = "739be429-bea8-5141-9913-cc70e7f3736d" -version = "1.0.3" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" - -[[deps.Measures]] -git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f" -uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" -version = "0.3.1" - -[[deps.Missings]] -deps = ["DataAPI"] -git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f" -uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.0.2" - -[[deps.Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" - -[[deps.NLSolversBase]] -deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "50310f934e55e5ca3912fb941dec199b49ca9b68" -uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.8.2" - -[[deps.NLsolve]] -deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] -git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1" -uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -version = "4.5.1" - -[[deps.NaNMath]] -git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe" -uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "0.3.7" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" - -[[deps.OffsetArrays]] -deps = ["Adapt"] -git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed" -uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.10.8" - -[[deps.Ogg_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" -uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" -version = "1.3.5+1" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" - -[[deps.OpenLibm_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "05823500-19ac-5b8b-9628-191a04bc5112" - -[[deps.OpenSSL_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "648107615c15d4e09f7eca16307bc821c1f718d8" -uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "1.1.13+0" - -[[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[deps.Optim]] -deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "bc0a748740e8bc5eeb9ea6031e6f050de1fc0ba2" -uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.6.2" - -[[deps.Opus_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720" -uuid = "91d4177d-7536-5919-b921-800302f37372" -version = "1.3.2+0" - -[[deps.OrderedCollections]] -git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.4.1" - -[[deps.PCRE_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488" -uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc" -version = "8.44.0+0" - -[[deps.Parameters]] -deps = ["OrderedCollections", "UnPack"] -git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" -uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" -version = "0.12.3" - -[[deps.Parsers]] -deps = ["Dates"] -git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa" -uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.2.2" - -[[deps.Pixman_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29" -uuid = "30392449-352a-5448-841d-b1acce4e97dc" -version = "0.40.1+0" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" - -[[deps.PlotThemes]] -deps = ["PlotUtils", "Requires", "Statistics"] -git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d" -uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" -version = "2.0.1" - -[[deps.PlotUtils]] -deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"] -git-tree-sha1 = "6f1b25e8ea06279b5689263cc538f51331d7ca17" -uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" -version = "1.1.3" - -[[deps.Plots]] -deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"] -git-tree-sha1 = "d16070abde61120e01b4f30f6f398496582301d6" -uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.25.12" - -[[deps.PositiveFactorizations]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" -uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" -version = "0.2.4" - -[[deps.Preferences]] -deps = ["TOML"] -git-tree-sha1 = "de893592a221142f3db370f48290e3a2ef39998f" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.4" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[deps.Qt5Base_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] -git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8" -uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1" -version = "5.15.3+0" - -[[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Random]] -deps = ["SHA", "Serialization"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[deps.Ratios]] -deps = ["Requires"] -git-tree-sha1 = "01d341f502250e81f6fec0afe662aa861392a3aa" -uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" -version = "0.4.2" - -[[deps.RealDot]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9" -uuid = "c1ae055f-0cd5-4b69-90a6-9a35b1a98df9" -version = "0.1.0" - -[[deps.RecipesBase]] -git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d" -uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -version = "1.2.1" - -[[deps.RecipesPipeline]] -deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"] -git-tree-sha1 = "995a812c6f7edea7527bb570f0ac39d0fb15663c" -uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" -version = "0.5.1" - -[[deps.Reexport]] -git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" -uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.2.2" - -[[deps.RelocatableFolders]] -deps = ["SHA", "Scratch"] -git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca" -uuid = "05181044-ff0b-4ac5-8273-598c1e38db00" -version = "0.1.3" - -[[deps.Requires]] -deps = ["UUIDs"] -git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" -uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.3.0" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" - -[[deps.Scratch]] -deps = ["Dates"] -git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda" -uuid = "6c6a2e73-6563-6170-7368-637461726353" -version = "1.1.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.SharedArrays]] -deps = ["Distributed", "Mmap", "Random", "Serialization"] -uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" - -[[deps.Showoff]] -deps = ["Dates", "Grisu"] -git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de" -uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" -version = "1.0.3" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[deps.SortingAlgorithms]] -deps = ["DataStructures"] -git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508" -uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.0.1" - -[[deps.SparseArrays]] -deps = ["LinearAlgebra", "Random"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - -[[deps.SpecialFunctions]] -deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.1.4" - -[[deps.Static]] -deps = ["IfElse"] -git-tree-sha1 = "65068e4b4d10f3c31aaae2e6cb92b6c6cedca610" -uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.5.6" - -[[deps.StaticArrays]] -deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "74fb527333e72ada2dd9ef77d98e4991fb185f04" -uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.4.1" - -[[deps.Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" - -[[deps.StatsAPI]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.2.1" - -[[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291" -uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.33.16" - -[[deps.StructArrays]] -deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"] -git-tree-sha1 = "57617b34fa34f91d536eb265df67c2d4519b8b98" -uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" -version = "0.6.5" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" - -[[deps.TableTraits]] -deps = ["IteratorInterfaceExtensions"] -git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" -uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" -version = "1.0.1" - -[[deps.Tables]] -deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"] -git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b" -uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.6.1" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[deps.URIs]] -git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355" -uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.3.0" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[deps.UnPack]] -git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b" -uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" -version = "1.0.2" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.UnicodeFun]] -deps = ["REPL"] -git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" -uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" -version = "0.4.1" - -[[deps.Unzip]] -git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670" -uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d" -version = "0.1.2" - -[[deps.Wayland_jll]] -deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23" -uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89" -version = "1.19.0+0" - -[[deps.Wayland_protocols_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da" -uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91" -version = "1.25.0+0" - -[[deps.WoodburyMatrices]] -deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3" -uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" -version = "0.5.5" - -[[deps.XML2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a" -uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.9.12+0" - -[[deps.XSLT_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] -git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" -uuid = "aed1982a-8fda-507f-9586-7b0439959a61" -version = "1.1.34+0" - -[[deps.Xorg_libX11_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] -git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527" -uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc" -version = "1.6.9+4" - -[[deps.Xorg_libXau_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e" -uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" -version = "1.0.9+4" - -[[deps.Xorg_libXcursor_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"] -git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd" -uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724" -version = "1.2.0+4" - -[[deps.Xorg_libXdmcp_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4" -uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05" -version = "1.1.3+4" - -[[deps.Xorg_libXext_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] -git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3" -uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3" -version = "1.3.4+4" - -[[deps.Xorg_libXfixes_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] -git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4" -uuid = "d091e8ba-531a-589c-9de9-94069b037ed8" -version = "5.0.3+4" - -[[deps.Xorg_libXi_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"] -git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246" -uuid = "a51aa0fd-4e3c-5386-b890-e753decda492" -version = "1.7.10+4" - -[[deps.Xorg_libXinerama_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"] -git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123" -uuid = "d1454406-59df-5ea1-beac-c340f2130bc3" -version = "1.1.4+4" - -[[deps.Xorg_libXrandr_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"] -git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631" -uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484" -version = "1.5.2+4" - -[[deps.Xorg_libXrender_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] -git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96" -uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa" -version = "0.9.10+4" - -[[deps.Xorg_libpthread_stubs_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb" -uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74" -version = "0.1.0+3" - -[[deps.Xorg_libxcb_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"] -git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6" -uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b" -version = "1.13.0+3" - -[[deps.Xorg_libxkbfile_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] -git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2" -uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a" -version = "1.1.0+4" - -[[deps.Xorg_xcb_util_image_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97" -uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b" -version = "0.4.0+1" - -[[deps.Xorg_xcb_util_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"] -git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1" -uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5" -version = "0.4.0+1" - -[[deps.Xorg_xcb_util_keysyms_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00" -uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7" -version = "0.4.0+1" - -[[deps.Xorg_xcb_util_renderutil_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e" -uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e" -version = "0.3.9+1" - -[[deps.Xorg_xcb_util_wm_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67" -uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361" -version = "0.4.1+1" - -[[deps.Xorg_xkbcomp_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"] -git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b" -uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4" -version = "1.4.2+4" - -[[deps.Xorg_xkeyboard_config_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"] -git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d" -uuid = "33bec58e-1273-512f-9401-5d533626f822" -version = "2.27.0+4" - -[[deps.Xorg_xtrans_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845" -uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" -version = "1.4.0+3" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" - -[[deps.Zstd_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927" -uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" -version = "1.5.2+0" - -[[deps.Zygote]] -deps = ["AbstractFFTs", "ChainRules", "ChainRulesCore", "DiffRules", "Distributed", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NaNMath", "Random", "Requires", "SparseArrays", "SpecialFunctions", "Statistics", "ZygoteRules"] -git-tree-sha1 = "93285d2877f1f1b09b2a2b029f90e9db10127022" -uuid = "e88e6eb3-aa80-5325-afca-941959d7151f" -version = "0.6.35" - -[[deps.ZygoteRules]] -deps = ["MacroTools"] -git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0" -uuid = "700de1a5-db45-46bc-99cf-38207098b444" -version = "0.2.2" - -[[deps.libass_jll]] -deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" -uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" -version = "0.15.1+0" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" - -[[deps.libfdk_aac_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" -uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" -version = "2.0.2+0" - -[[deps.libpng_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c" -uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" -version = "1.6.38+0" - -[[deps.libvorbis_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] -git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c" -uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" -version = "1.3.7+1" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" - -[[deps.x264_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" -uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" -version = "2021.5.5+0" - -[[deps.x265_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" -uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" -version = "3.5.0+0" - -[[deps.xkbcommon_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"] -git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6" -uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd" -version = "0.9.1+5" diff --git a/README.md b/README.md index c260935..696ba7c 100644 --- a/README.md +++ b/README.md @@ -6,18 +6,20 @@ This package implements a variety of models for rates of electrochemical reactio ### Building and Evaluating Models We can construct a simple Butler-Volmer model with symmetric electron transfer as follows: ```julia -julia> bv = ButlerVolmer(1.0) +julia> bv = ButlerVolmer() ButlerVolmer(A=1.0, α=0.5) ``` To evaluate the model and find the reaction rate at a given voltage is as simple as: ```julia julia> bv(0.1) -6.695821798442868 +6.865409234059785 ``` +(note that this is equivalent to `rate_constant(0.1, bv)`, but this "callable" syntax is often convenient) + By default, this computes the net reaction rate. We can also get the rate only in the oxidative direction by passing the `ox` flag: ```julia julia> bv(0.1, true) -6.84197835551441 +7.0081012366762145 ``` Let's try a fancier model now. The `MarcusHushChidseyDOS` model requires information about the density of states of the electrode, which we can read from a file: ```julia @@ -26,7 +28,7 @@ MarcusHushChidseyDOS(A=20.0, λ=0.2) ``` In this case, since `MarcusHushChidseyDOS<:IntegralModel`, we can't directly evaluate the object as a callable because it's ambiguous if we want the integrand or the full integrated rate constant, and we must specify: ```julia -julia> compute_k(0.1, mhcd) +julia> rate_constant(0.1, mhcd) 3.426690873746331 ``` By default, the integration bounds are chosen to be the energy bounds over which the DOS is defined, which we can easily check: @@ -42,8 +44,8 @@ Another feature we support for `MarcusHushChidseyDOS` models is incorporating th julia> mhcd = MarcusHushChidseyDOS(0.82, "data/DOSes/bernal_graphene.txt") MarcusHushChidseyDOS(A=1.0, λ=0.82) -julia> compute_k(0.4, mhcd; calc_cq=true, Vq_min=-0.45, Vq_max=0.45) -0.00044958691359217384 +julia> rate_constant(0.15, mhcd, true; calc_cq=true, Vq_min=-0.3, Vq_max=0.3, E_min = -0.6, E_max=0.6) +2.8017566534572208e-5 ``` ### Fitting and Comparing Models @@ -77,43 +79,44 @@ julia> plot_models(models) ``` models +### Electrochemical Phase Diagrams +For more on this, see upcoming paper: [arxiv link placeholder] + ## Supported Models -(I hope the equations below are visible in both dark and light mode...) ### Butler-Volmer -Probably the most basic (and entirely empirical) kinetic model. The rate constants are given by: +Probably the most basic (and largely empirical) kinetic model. The rate constants are given by: -BV_ox +$$k_{\text{ox}}^{\text{BV}}(\eta)=A\exp\left(\frac{\alpha\eta}{k_{\text B}T}\right)\\k_{\text{red}}^{\text{BV}}(\eta)=A\exp\left(\frac{(1-\alpha)\eta}{k_{\text B}T}\right)$$ -BV_red - -Implemented in the package as `ButlerVolmer`, possessing two fields: the prefactor `A` and the transfer coefficient `α`. +It is implemented in the package as `ButlerVolmer`, possessing two fields: the prefactor `A` and the transfer coefficient `α`. ### Marcus One step up with a bit of mechanism... -Marcus_k +$$k_{\text{red/ox}}^{\text{Marcus}}(\eta) = A\exp\left(-\frac{(\lambda\pm\eta)^2}{4k_{\text B}T}\right)$$ Implemented in the package as `Marcus`, with parameters for prefactor `A` and reorganization energy `λ`. ### Marcus-Hush-Chidsey Originated from [Chidsey's 1991 paper](https://dx.doi.org/10.1126/science.251.4996.919). -MHC_k +$$k_{\text{red/ox}}^{\text{MHC}}(\eta) = A\int_{-\infty}^\infty \exp\left(-\frac{(\varepsilon-\lambda\mp\eta)^2}{4k_{\text B}T}\right)f_{\text{FD}}(\varepsilon, T)\mathrm d\varepsilon$$ Implemented as `MarcusHushChidsey`, with prefactor `A`, reorganization energy `λ`, and also an `average_dos` parameter (which defaults to 1.0), which is useful for making direct comparisons against `MarcusHushChidseyDOS` (see below). ### Asymptotic Marcus-Hush-Chidsey -Originated from [Zeng et al. 2014](https://dx.doi.org/10.1016/j.jelechem.2014.09.038d) (equation 17). Note that in their notation, $\eta$ is scaled by the thermal energy. +Originated from [Zeng et al. 2014](https://dx.doi.org/10.1016/j.jelechem.2014.09.038d) (equation 17). We have corrected it to recover the correct temperature scaling, [paper link placeholder] Implemented as `AsymptoticMarcusHushChidsey` with parameters `A` and `λ`. ### Marcus-Hush-Chidsey + DOS Originated from [Kurchin and Viswanathan 2020](https://dx.doi.org/10.1063/5.0023611 ). -MHCD_ox -MHCD_red +$$k_{\text{ox}}^{\text{MHCKV}}(\eta)=A\int_{-\infty}^\infty\mathcal{D}(\varepsilon)\exp\left(-\frac{(\lambda-\eta+\varepsilon)^2}{4\lambda k_{\text{B}}T}\right)(1-f_{\text{FD}}(\varepsilon,T))\mathrm{d}\varepsilon\\ +k_{\text{red}}^{\text{MHCKV}}(\eta)=A\int_{-\infty}^\infty\mathcal{D}(\varepsilon)\exp\left(-\frac{(\lambda+\eta-\varepsilon)^2}{4\lambda k_{\text{B}}T}\right)f_{\text{FD}}(\varepsilon,T)\mathrm{d}\varepsilon$$ + -Implemented as `MarcusHushChidseyDOS`, with parameters `A`, `λ`, and `dos`, which must be a `DOSData` object. +Implemented as `MarcusHushChidseyDOS`, with parameters `A`, `λ`, and `dos`, which must be a `DOSData` object (can be automatically constructed from a two-column text file of energies and DOS values). ## Contributing We are happy to take pull requests for new features, new models, etc. by pull request! It is suggested (though not required) that you first open an issue to discuss. diff --git a/benchmark/bench_suite.jl b/benchmark/bench_suite.jl index cfd566b..46604f0 100644 --- a/benchmark/bench_suite.jl +++ b/benchmark/bench_suite.jl @@ -8,7 +8,7 @@ suite["integrals"] = BenchmarkGroup() suite["integrals"]["scale"] = @benchmarkable ElectrochemicalKinetics.scale(0.1, 0.2) seconds=time_mult # reusing the same ones from phase_diagram_tests -bv = ButlerVolmer(300) +bv = ButlerVolmer(300, 0.5) m = Marcus(5000, 0.3) amhc = AsymptoticMarcusHushChidsey(70000, 0.3) ni_kms = [bv, m, amhc] @@ -22,8 +22,8 @@ suite["rate constants"]["non-integral models"]["vector V"] = BenchmarkGroup() for km in ni_kms modeltype = typeof(km) - suite["rate constants"]["non-integral models"]["scalar V"][modeltype] = @benchmarkable compute_k(0.1, $km) seconds=time_mult*0.1 - suite["rate constants"]["non-integral models"]["vector V"][modeltype] = @benchmarkable compute_k(0.01:0.01:0.2, $km) seconds=time_mult + suite["rate constants"]["non-integral models"]["scalar V"][modeltype] = @benchmarkable rate_constant($0.1, $km) seconds=time_mult*0.1 + suite["rate constants"]["non-integral models"]["vector V"][modeltype] = @benchmarkable rate_constant($0.01:0.01:0.2, $km) seconds=time_mult end mhc = MarcusHushChidsey(70000, 0.3, 1.0) @@ -38,11 +38,11 @@ suite["rate constants"]["integral models"]["vector V"] = BenchmarkGroup() for km in i_kms modeltype = names[km] - suite["rate constants"]["integral models"]["scalar V"][modeltype] = @benchmarkable compute_k(0.1, $km) seconds=time_mult*8 - suite["rate constants"]["integral models"]["vector V"][modeltype] = @benchmarkable compute_k(0.01:0.01:0.2, $km) seconds=time_mult*11 + suite["rate constants"]["integral models"]["scalar V"][modeltype] = @benchmarkable rate_constant($0.1, $km) seconds=time_mult*8 + suite["rate constants"]["integral models"]["vector V"][modeltype] = @benchmarkable rate_constant($0.01:0.01:0.2, $km) seconds=time_mult*11 end -# fit_overpotential benchmarks +# overpotential benchmarks # TODO: add vector cases k = 2e3 @@ -53,8 +53,8 @@ suite["overpotential fitting"]["non-integral models"]["AD"] = BenchmarkGroup() for km in ni_kms modeltype = typeof(km) - suite["overpotential fitting"]["non-integral models"]["non-AD"][modeltype] = @benchmarkable fit_overpotential($km, $k, autodiff=false) seconds=time_mult*3 - suite["overpotential fitting"]["non-integral models"]["AD"][modeltype] = @benchmarkable fit_overpotential($km, $k, autodiff=true) seconds=time_mult*5 + suite["overpotential fitting"]["non-integral models"]["non-AD"][modeltype] = @benchmarkable overpotential($k, $km, autodiff=false) seconds=time_mult*3 + suite["overpotential fitting"]["non-integral models"]["AD"][modeltype] = @benchmarkable overpotential($k, $km, autodiff=true) seconds=time_mult*5 end suite["overpotential fitting"]["integral models"] = BenchmarkGroup() @@ -62,6 +62,16 @@ suite["overpotential fitting"]["integral models"]["non-AD"] = BenchmarkGroup() suite["overpotential fitting"]["integral models"]["AD"] = BenchmarkGroup() for km in i_kms[2:end] # skipping MHC for now modeltype = names[km] - suite["overpotential fitting"]["integral models"]["non-AD"][modeltype] = @benchmarkable fit_overpotential($km, $k, autodiff=false) seconds=time_mult*10 - suite["overpotential fitting"]["integral models"]["AD"][modeltype] = @benchmarkable fit_overpotential($km, $k, autodiff=true) seconds=time_mult*20 + suite["overpotential fitting"]["integral models"]["non-AD"][modeltype] = @benchmarkable overpotential($k, $km, autodiff=false) seconds=time_mult*10 + suite["overpotential fitting"]["integral models"]["AD"][modeltype] = @benchmarkable overpotential($k, $km, autodiff=true) seconds=time_mult*20 end + +suite["phase diagrams"] = BenchmarkGroup() +suite["phase diagrams"]["non-integral models"] = BenchmarkGroup() +for km in ni_kms + modeltype = typeof(km) + suite["phase diagrams"]["non-integral models"][modeltype] = @benchmarkable phase_diagram($km, T=350) + # TODO: add a lower-temp BV one +end + +# TODO: get working, add benchmarks for integral model phase diagrams diff --git a/benchmark/results_branch.jls b/benchmark/results_branch.jls new file mode 100644 index 0000000..b5e7933 Binary files /dev/null and b/benchmark/results_branch.jls differ diff --git a/benchmark/results_master.jls b/benchmark/results_master.jls new file mode 100644 index 0000000..1d87362 Binary files /dev/null and b/benchmark/results_master.jls differ diff --git a/img/MHC.png b/img/MHC.png deleted file mode 100644 index 0af815e..0000000 Binary files a/img/MHC.png and /dev/null differ diff --git a/img/MHCD_ox.png b/img/MHCD_ox.png deleted file mode 100644 index 2db79f6..0000000 Binary files a/img/MHCD_ox.png and /dev/null differ diff --git a/img/MHCD_red.png b/img/MHCD_red.png deleted file mode 100644 index adc8c58..0000000 Binary files a/img/MHCD_red.png and /dev/null differ diff --git a/img/Marcus.png b/img/Marcus.png deleted file mode 100644 index 565e873..0000000 Binary files a/img/Marcus.png and /dev/null differ diff --git a/img/bv_ox.png b/img/bv_ox.png deleted file mode 100644 index 240d4be..0000000 Binary files a/img/bv_ox.png and /dev/null differ diff --git a/img/bv_red.png b/img/bv_red.png deleted file mode 100644 index 75539d4..0000000 Binary files a/img/bv_red.png and /dev/null differ diff --git a/img/plot1.png b/img/plot1.png index ebe5318..a231c76 100644 Binary files a/img/plot1.png and b/img/plot1.png differ diff --git a/img/plot2.png b/img/plot2.png index ff4c14c..49d96bf 100644 Binary files a/img/plot2.png and b/img/plot2.png differ diff --git a/src/ElectrochemicalKinetics.jl b/src/ElectrochemicalKinetics.jl index a465a17..52fe134 100644 --- a/src/ElectrochemicalKinetics.jl +++ b/src/ElectrochemicalKinetics.jl @@ -1,29 +1,28 @@ module ElectrochemicalKinetics -include("dos.jl") +const kB = 8.61733326e-5 + +include("utils/dos.jl") export DOS using .DOS: DOSData, get_dos export DOSData, get_dos -include("kinetic_models.jl") +include("kinetic_models/kinetic_models.jl") export ButlerVolmer, AsymptoticMarcusHushChidsey export fermi_dirac export integrand -export KineticModel, IntegralModel +export KineticModel, NonIntegralModel, IntegralModel export Marcus, MarcusHushChidsey, MarcusHushChidseyDOS - -include("quantum_capacitance.jl") - -include("rate_constant.jl") -export compute_k, compute_k_cq +export rate_constant, rate_constant_cq include("fitting.jl") -include("integration_utils.jl") -export fitting_params, fit_model, fit_overpotential, is_dosmodel +include("utils/integration_utils.jl") +export fitting_params, fit_model, overpotential, is_dosmodel include("phase_diagrams.jl") export kB, h, s, g_thermo, µ_thermo -export µ_kinetic, g_kinetic, common_tangent, find_phase_boundaries +export µ_kinetic, g_kinetic +export common_tangent, find_phase_boundaries, phase_diagram include("plot_fcns.jl") export plot_models, plot_exp_and_models diff --git a/src/fitting.jl b/src/fitting.jl index f843eaa..64fff02 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -9,75 +9,126 @@ log_loss(y, y_pred) = (log.(y) .- log.(y_pred)).^2 # sum of squares loss in linear coordinates linear_loss(y, y_pred) = (y .- y_pred).^2 +# one that works for both signs +function janky_log_loss(y, y_pred) + if all(y .* y_pred .> 0) + return log_loss(abs.(y), abs.(y_pred)) + else + # @warn "Values have different signs...defaulting to squared error" + # this should be fine to (hopefully) push things into the right region... + return linear_loss(y, y_pred) + end +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))) + # make sure guess has correct sign + if k isa Real + if k == 0 + guess .= 0 + elseif k < 0 + guess[guess .> 0] .= -1 * guess[guess .> 0] + end + elseif k isa Vector + guess[k .== 0] .= 0 + guess[k .< 0 .&& guess .> 0] .= -1 * guess[k .< 0 .&& guess .> 0] + end + return guess +end + +# TODO: support reaction direction flag here """ - fit_overpotential(model, k; kwargs...) + overpotential(k, model; kwargs...) -Given values for current/rate constant and specified model parameters, find the overpotential that would have resulted in it. (This is the inverse of the `compute_k` function.) +Given values for current/rate constant and specified model parameters, find the overpotential that would have resulted in it. (This is the inverse of the `rate_constant` function.) + +NOTE that this currently only solves for net reaction rates. """ -function fit_overpotential(model::KineticModel, k, guess = fill(0.1, length(k)); kT = 0.026, loss = log_loss, autodiff = true, verbose=false, kwargs...) - function compare_k!(storage, V) - storage .= loss(k, compute_k(V, model; kT = kT, kwargs...)) +function overpotential(k, model::KineticModel, guess = _get_guess(k, model); T = 298, loss = janky_log_loss, autodiff = true, verbose=false, kwargs...) + # wherever k=0 we can shortcut since the answer has to be 0 + k_solve = k + if k==0 # scalar k=0, possibly vector model + if verbose + println("shortcutting full solve") + end + return zeros(Float32, length(model)) + elseif any(k .== 0) # vector k with some zero elements, scalar model + if verbose + println("shortcutting part of solve") + end + k_solve = k[k.!=0] end + function compare_k!(storage, V) + storage .= loss(k, rate_constant(V, model; T = T, kwargs...)) + end Vs = if autodiff function myfun!(F, J, V) - if J == nothing - F .= loss(k, compute_k(V, model; kT = kT, kwargs...)) - elseif F == nothing && !isnothing(J) - # TODO: we could speed up the case of scalar k's by dispatching to call gradient here instead - gs = Zygote.gradient(V[1]) do 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...)) + elseif isnothing(F) && !isnothing(J) + gs = Zygote.gradient(V) do V Zygote.forwarddiff(V) do V - loss(k, compute_k(V, model; kT = kT, kwargs...)) |> sum + loss(k_solve, rate_constant(V, model; T = T, kwargs...)) |> sum end end[1] - J .= gs + J .= diagm(gs) else - y, back = Zygote.pullback(V[1]) do V + y, back = Zygote.pullback(V) do V Zygote.forwarddiff(V) do V - loss(k, compute_k(V, model; kT = kT, kwargs...)) |> sum + loss(k_solve, rate_constant(V, model; T = T, kwargs...)) |> sum end end F .= y - J .= back(one.(y))[1] + J .= diagm(back(one.(y))[1]) end + # println("J=", J) end - Vs = nlsolve(only_fj!(myfun!), guess, show_trace=verbose) + Vs = nlsolve(only_fj!(myfun!), guess, show_trace=verbose, extended_trace=true) else Vs = nlsolve(compare_k!, guess, show_trace=verbose) end if !converged(Vs) @warn "Overpotential fit not fully converged...you may have fed in an unreachable reaction rate!" end - # this is janky but the broadcast doesn't work otherwise - if typeof(k) <: Array - Vs.zero + sol = Vs.zero + sol_return = zero(guess) + if k_solve != k # need to put back the zero elements + sol_return[k.!=0] .= sol else - Vs.zero[1] + sol_return = sol 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) -Zygote.@adjoint function fit_overpotential(model, k, guess; loss = log_loss, kT = 0.026, autodiff = true, verbose = false, kw...) - Vs = fit_overpotential(model, k, guess; loss, verbose, autodiff, kT, kw...) +# TODO: test this +Zygote.@adjoint function overpotential(k, model, guess; loss = log_loss, T = 298, autodiff = true, verbose = false, kw...) + Vs = overpotential(k, model, guess; loss=loss, verbose=verbose, autodiff=autodiff, T=T, kw...) function back(vs) gs = Zygote.jacobian(vs) do V Zygote.forwarddiff(V) do V - compute_k(V, model; kw...) + rate_constant(V, model; T=T, kw...) end end[1] inv.(gs) end - Vs, Δ -> (nothing, nothing, Δ .* back(Vs)) + Vs, Δ -> (Δ .* back(Vs), nothing, nothing) end -# basically the gradient of fit_overpotential should just be the inverse of the gradient of the forward solve at that point, i.e. if compute_k(V, model) ==k then fit_overpotential(model, k)==V , so we don’t need to diff through the actual solve in fit_overpotential because the gradient should just be the inverse of the gradient of compute_k at V +# basically the gradient of overpotential should just be the inverse of the gradient of the forward solve at that point, i.e. if rate_constant(V, model) ==k then overpotential(model, k)==V , so we don’t need to diff through the actual solve in overpotential because the gradient should just be the inverse of the gradient of rate_constant at V # multiple models, one k value (used in thermo example) -fit_overpotential(models::Vector{<:KineticModel}, k::Real, guess=fill(0.1, length(k)); kwargs...) = fit_overpotential.(models, Ref(k), 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{<:KineticModel}) = fieldnames(t) +fitting_params(t::Type{<:NonIntegralModel}) = fieldnames(t) fitting_params(::Type{MarcusHushChidsey}) = (:A, :λ) fitting_params(::Type{MarcusHushChidseyDOS}) = (:A, :λ) const default_param_bounds = Dict(:A => (0.1, 50000), :λ => (0.01, 0.5), :α => (0.01, 0.99)) @@ -92,7 +143,7 @@ const default_param_bounds = Dict(:A => (0.1, 50000), :λ => (0.01, 0.5), :α => # Keyword Arguments Requirements differ by model type... * `ButlerVolmer`, `AsymptoticMarcusHushChidsey`, `Marcus`: none -* `MarcusHushChidsey`: average_dos::Float64 OR dos::DOSData OR dos_file::String +* `MarcusHushChidsey`: average_dos::Float32 OR dos::DOSData OR dos_file::String * `MarcusHushChidseyDOS`: dos::DOSData OR dos_file Some are always options... * `param_bounds::Dict{Symbol,Any}`: ranges of guesses for relevant model parameters. (must include all necessary keys, but defaults to some sensible ranges if not provided, see `default_param_bounds`...note that you should provide this for faster fitting if you know bounds) @@ -102,11 +153,11 @@ function fit_model( exp_data::Matrix, model_type::Type{<:KineticModel}; param_bounds::Dict = default_param_bounds, - kT::Real = 0.026, + T = 298, kwargs... ) V_vals = exp_data[:, 1] - eval_model(model) = [compute_k(V, model; kT = kT) for V in V_vals] + eval_model(model) = [rate_constant(V, model; T=T) for V in V_vals] _fit_model(exp_data, model_type, param_bounds, eval_model) end @@ -115,15 +166,15 @@ function fit_model( exp_data::Matrix, model_type::Type{<:IntegralModel}; param_bounds::Dict = default_param_bounds, - kT::Real = 0.026, - E_min = -100 * kT, - E_max = 100 * kT, + T = 298, + E_min = -100 * kB*T, + E_max = 100 * kB*T, kwargs..., ) V_vals = exp_data[:, 1] - eval_model(model) = compute_k(V_vals, + eval_model(model) = rate_constant(V_vals, model; - kT = kT, + T = T, E_min = E_min, E_max = E_max, kwargs...) @@ -140,7 +191,7 @@ function _fit_model( I_vals = exp_data[:, 2] model_builder = _get_model_builder(model_type, param_bounds; kwargs...) - sq_error(I_pred) = sum((I_pred .- I_vals) .^ 2) + sq_error(I_pred) = sum((abs.(I_pred) .- I_vals) .^ 2) # WTF, why do I need these here again...just when I think I understand scope fitting_params(t::Type{<:KineticModel}) = fieldnames(t) @@ -148,7 +199,6 @@ function _fit_model( fitting_params(::Type{MarcusHushChidseyDOS}) = (:A, :λ) # find best-fitting params - # Zygote is tripped up by QuadGK, so that one has to be done with black-box optimization, but the non-integral models work with autodiff opt_func = params -> sq_error(model_evaluator(model_builder(params))) local best_params # function grad!(s, x) @@ -157,18 +207,18 @@ function _fit_model( # s[i] = gs[i] # end # end - lower = Float64.([param_bounds[p][1] for p in fitting_params(model_type)]) - upper = Float64.([param_bounds[p][2] for p in fitting_params(model_type)]) + lower = Float32.([param_bounds[p][1] for p in fitting_params(model_type)]) + upper = Float32.([param_bounds[p][2] for p in fitting_params(model_type)]) init_guess = 0.5 .* (lower .+ upper) # set optimisers based on the model type # IntegralModels performed better with GradientDescent optimizer, opts = if model_type <: IntegralModel opts = Optim.Options(show_trace = false, - iterations = 3, - outer_iterations = 4, + iterations = 5, + outer_iterations = 7, x_tol = 1., - f_tol = 1e3) + f_tol = 1e2) Fminbox(NelderMead()), opts else Fminbox(NelderMead()), Optim.Options() diff --git a/src/kinetic_models.jl b/src/kinetic_models.jl deleted file mode 100644 index 20c7525..0000000 --- a/src/kinetic_models.jl +++ /dev/null @@ -1,240 +0,0 @@ -using .DOS -using SpecialFunctions - -""" - fermi_dirac(E, kT=0.026) - -Compute the value of the Fermi-Dirac distribution at energy `E` (relative to the Fermi energy) and thermal energy `kT`. -""" -fermi_dirac(E; kT = 0.026) = inv.(1 .+ exp.(E ./ kT)) - -abstract type KineticModel end - -# dispatch for net rates, returns absolute value -(km::KineticModel)(V_app; kT = 0.026) = - abs.(km(V_app, Val(true); kT = kT) - km(V_app, Val(false); kT = kT)) - -# return a new one with a scaled prefactor -import Base.* -function *(c::Real, km::KineticModel) - new_A = c*km.A - other_params = getfield.([km], propertynames(km))[2:end] - typeof(km)(new_A, other_params...) -end - -# generic pretty printing -function Base.show(io::IO, m::KineticModel) - s = repr(typeof(m)) * "(" - for field in propertynames(m) - s = - s * - string(field) * - "=" * - string(round(getproperty(m, field), sigdigits = 3)) * - ", " - end - s = s[1:end-2] * ")" - print(io, s) -end - -""" - ButlerVolmer(A, α) - ButlerVolmer(A) - -Computes Butler-Volmer kinetics. - -If initialized with one argument, assumes symmetric electron transfer (α=0.5) and sets this to be the prefactor A. Note that this prefactor implicitly contains information about equilibrium activation energies, as well as geometric information. -""" -struct ButlerVolmer <: KineticModel - A::Float64 - α::Float64 -end - -# default to unit prefactor and symmetric response -ButlerVolmer() = ButlerVolmer(1.0, 0.5) -ButlerVolmer(A) = ButlerVolmer(A, 0.5) - -function (bv::ButlerVolmer)(V_app, ::Val{true}; kT::Real = 0.026) - exp_arg = (bv.α .* V_app) ./ kT - bv.A .* exp.(exp_arg) -end - -function (bv::ButlerVolmer)(V_app, ::Val{false}; kT::Real = 0.026) - exp_arg = -((1 - bv.α) .* V_app) ./ kT - bv.A .* exp.(exp_arg) -end - -""" - Marcus(A, λ) - Marcus(λ) - -Computes Marcus kinetics. - -If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor A to 1.0. -""" -struct Marcus <: KineticModel - A::Float64 - λ::Float64 -end - -# default prefactor is 1 -Marcus(λ) = Marcus(1.0, λ) - -function (m::Marcus)(V_app, ox::Bool; kT::Real = 0.026) - local exp_arg - if ox - exp_arg = -(m.λ .+ V_app).^2 ./ (4 * m.λ * kT) - else - exp_arg = -(m.λ .- V_app).^2 ./ (4 * m.λ * kT) - end - m.A .* exp.(exp_arg) -end - -function (m::Marcus)(V_app, ::Val{true}; kT::Real = 0.026) - exp_arg = -(m.λ .+ V_app).^2 ./ (4 * m.λ * kT) - m.A .* exp.(exp_arg) -end - -function (m::Marcus)(V_app, ::Val{false}; kT::Real = 0.026) - exp_arg = -(m.λ .- V_app).^2 ./ (4 * m.λ * kT) - m.A .* exp.(exp_arg) -end - -""" - AsymptoticMarcusHushChidsey(A, λ) - AsymptoticMarcusHushChidsey(λ) - -Computes asymptotic solution to MHC model, as described in Zeng et al.: 10.1016/j.jelechem.2014.09.038d, with a corrction prefactor of kT since there is an error in the nondimensionalization in that work. - -If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor to 1.0. -""" -struct AsymptoticMarcusHushChidsey <: KineticModel - A::Float64 - λ::Float64 -end - -# default prefactor is 1 -AsymptoticMarcusHushChidsey(λ) = AsymptoticMarcusHushChidsey(1.0, λ) - -function (amhc::AsymptoticMarcusHushChidsey)(V_app, ox::Bool; kT::Real = 0.026) - η = (2 * ox - 1) .* V_app ./ kT - λ_nondim = amhc.λ / kT - a = 1 + sqrt(λ_nondim) - arg = (λ_nondim .- sqrt.(a .+ η.^2)) ./ (2 * sqrt(λ_nondim)) - pref = sqrt(π * λ_nondim) ./ (1 .+ exp.(-η)) - return kT * amhc.A .* pref .* erfc.(arg) -end - - -function (amhc::AsymptoticMarcusHushChidsey)(V_app; kT::Real = 0.026) - η = V_app ./ kT - λ_nondim = amhc.λ / kT - a = 1 + sqrt(λ_nondim) - arg = (λ_nondim .- sqrt.(a .+ η.^2)) ./ (2 * sqrt(λ_nondim)) - pref = sqrt(π * λ_nondim) .* tanh.(η/2) - return abs.(kT * amhc.A .* pref .* erfc.(arg)) -end - - -abstract type IntegralModel <: KineticModel end # "Marcus-like" - -# check to catch missed dispatches for new types -# 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 -# dispatch for net rates -integrand(km::IntegralModel, V_dl; kwargs...) = - E -> abs.( - integrand(km, V_dl, true; kwargs...)(E) - integrand(km, V_dl, false; kwargs...)(E), - ) - -""" - MarcusHushChidsey(A, λ, average_dos) - MarcusHushChidsey(λ, average_dos) - MarcusHushChidsey(λ) - -Computes Marcus-Hush-Chidsey kinetics: 10.1126/science.251.4996.919 - -Note that for "typical" reorganization energy values (in the vicinity of 10*kT at typical temperatures, i.e. a few tenths of an eV), `AsymptoticMarcusHushChidsey` is comparably accurate to and much faster to evaluate than this model. - -If either the prefactor or the average dos are omitted, their values are assumed to be 1. Note that strictly speaking, `average_dos` and the prefactor `A` are redundant. They are both included primarily to facilitate comparisons with similarly parametrized Marcus-like models such as `MarcusHushChidseyDOS`. -""" -struct MarcusHushChidsey <: IntegralModel - A::Real - λ::Real - average_dos::Real -end - -# default prefactor is 1 -MarcusHushChidsey(λ, avg_dos) = MarcusHushChidsey(1.0, λ, avg_dos) -# assume prefactor = 1 and avg_dos = 1 -MarcusHushChidsey(λ) = MarcusHushChidsey(1.0, λ, 1.0) -# convert more detailed DOS information to just pull out average -MarcusHushChidsey(A, λ, dd::DOSData) = MarcusHushChidsey(A, λ, dd.average_value) -MarcusHushChidsey(A, λ, dos_file::Union{Matrix,String}; kwargs...) = - MarcusHushChidsey(A, λ, DOSData(dos_file; kwargs...)) - -# TODO: Check that both this and +DOS versions still match original paper -function integrand(mhc::MarcusHushChidsey, V_dl, ox::Bool; kT::Real = 0.026) - function marcus_term(E) - local exp_arg - if ox - exp_arg = -((E .- mhc.λ .+ V_dl) .^ 2) ./ (4 * mhc.λ * kT) - else - exp_arg = -((E .- mhc.λ .- V_dl) .^2) ./ (4 * mhc.λ * kT) - end - exp.(exp_arg) - end - f(E::Vector) = hcat((mhc.A * mhc.average_dos) .* marcus_term.(E) .* fermi_dirac.(E; kT = kT)...)' # will return a matrix of both E and V_dl are vectors, first index will be energies and second voltages - f(E::Real) = (mhc.A * mhc.average_dos) .* marcus_term.(E) .* fermi_dirac.(E; kT = kT) -end - - - -""" - MarcusHushChidseyDOS(A=1.0, λ, dos) - MarcusHushChidseyDOS(A=1.0, λ, dos_file) - -Computes Marcus-Hush-Chidsey + DOS kinetics as described in Kurchin and Viswanathan: 10.1063/5.0023611 -""" -struct MarcusHushChidseyDOS <: IntegralModel - A::Float64 - λ::Float64 - dos::DOSData -end - -function Base.show(io::IO, mhcd::MarcusHushChidseyDOS) - s = repr(typeof(mhcd)) * "(" - s *= "A=$(round(mhcd.A, sigdigits=3)), λ=$(round(mhcd.λ, sigdigits=3)))" - print(io, s) -end - -# default prefactor is 1 -MarcusHushChidseyDOS(λ, dd::DOSData) = MarcusHushChidseyDOS(1.0, λ, dd) - -MarcusHushChidseyDOS(A, λ, dos_file::Union{Matrix,String}; kwargs...) = - MarcusHushChidseyDOS(A, λ, DOSData(dos_file; kwargs...)) - -MarcusHushChidseyDOS(λ, dos_file::Union{Matrix,String}; kwargs...) = MarcusHushChidseyDOS(1.0, λ, DOSData(dos_file; kwargs...)) - -function integrand( - mhcd::MarcusHushChidseyDOS, - V_dl, - ox::Bool; - kT::Real = 0.026, - V_q = 0.0, -) - function marcus_term(E) - local exp_arg - if ox - exp_arg = -(( (mhcd.λ .- V_dl) .+ E) .^ 2) ./ (4 * mhcd.λ * kT) - else - exp_arg = -(( (mhcd.λ .+ V_dl) .- E) .^ 2) ./ (4 * mhcd.λ * kT) - end - exp.(exp_arg) - end - fd(E) = ox ? 1 .- fermi_dirac(E; kT = kT) : fermi_dirac(E; kT = kT) - # TODO: add two dispatches as with MHC above - E -> mhcd.A .* mhcd.dos.interp_func.(E .+ V_q) .* fd(E) .* marcus_term(E) -end diff --git a/src/kinetic_models/AsymptoticMarcusHushChidsey.jl b/src/kinetic_models/AsymptoticMarcusHushChidsey.jl new file mode 100644 index 0000000..0ec51c6 --- /dev/null +++ b/src/kinetic_models/AsymptoticMarcusHushChidsey.jl @@ -0,0 +1,38 @@ +""" + AsymptoticMarcusHushChidsey(A, λ) + AsymptoticMarcusHushChidsey(λ) + +Computes asymptotic solution to MHC model, as described in Zeng et al.: 10.1016/j.jelechem.2014.09.038d, with a corrction prefactor of kT since there is an error in the nondimensionalization in that work. + +If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor to 1.0. +""" +struct AsymptoticMarcusHushChidsey{T} <: NonIntegralModel{T} + A::T + λ::T + function AsymptoticMarcusHushChidsey(A, λ) + ps = consistent_params(Float32.(A), Float32.(λ)) + new{typeof(ps[1])}(ps...) + end +end + +# default prefactor is 1 +AsymptoticMarcusHushChidsey(λ) = AsymptoticMarcusHushChidsey(1.0, λ) + +function rate_constant(V_app, amhc::AsymptoticMarcusHushChidsey, ox::Bool; T = 298) + η = (2 * ox - 1) .* V_app ./ (kB * T) + λ_nondim = amhc.λ / (kB * T) + a = 1 .+ sqrt.(λ_nondim) + arg = (λ_nondim .- sqrt.(a .+ η.^2)) ./ (2 .* sqrt.(λ_nondim)) + pref = sqrt.(π .* λ_nondim) ./ (1 .+ exp.(-η)) + 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 diff --git a/src/kinetic_models/ButlerVolmer.jl b/src/kinetic_models/ButlerVolmer.jl new file mode 100644 index 0000000..27539ab --- /dev/null +++ b/src/kinetic_models/ButlerVolmer.jl @@ -0,0 +1,31 @@ +""" + ButlerVolmer(A, α) + ButlerVolmer(A) + +Computes Butler-Volmer kinetics. + +If initialized with one argument, assumes symmetric electron transfer (α=0.5) and sets this to be the prefactor A. Note that this prefactor implicitly contains information about equilibrium activation energies, as well as geometric information. +""" +struct ButlerVolmer{T} <: NonIntegralModel{T} + A::T + α::T + function ButlerVolmer(A, α) + @assert all(α .<= 1.0 .&& α .>=0.0) "Electron transfer coefficient must be in [0,1]" + ps = consistent_params(Float32.(A), Float32.(α)) + new{typeof(ps[1])}(ps...) + end +end + +# default to unit prefactor and symmetric response +ButlerVolmer() = ButlerVolmer(1.0, 0.5) +ButlerVolmer(α) = ButlerVolmer(1.0, α) + +function rate_constant(V_app, bv::ButlerVolmer, ::Val{true}; T = 298) + exp_arg = (bv.α .* V_app) / (kB * T) + bv.A .* exp.(exp_arg) +end + +function rate_constant(V_app, bv::ButlerVolmer, ::Val{false}; T = 298) + exp_arg = -((1 .- bv.α) .* V_app) / (kB * T) + bv.A .* exp.(exp_arg) +end diff --git a/src/kinetic_models/Marcus.jl b/src/kinetic_models/Marcus.jl new file mode 100644 index 0000000..3697df1 --- /dev/null +++ b/src/kinetic_models/Marcus.jl @@ -0,0 +1,29 @@ +""" + Marcus(A, λ) + Marcus(λ) + +Computes Marcus kinetics. + +If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor A to 1.0. +""" +struct Marcus{T} <: NonIntegralModel{T} + A::T + λ::T + function Marcus(A, λ) + ps = consistent_params(Float32.(A), Float32.(λ)) + new{typeof(ps[1])}(ps...) + end +end + +# default prefactor is 1 +Marcus(λ) = Marcus(1.0, λ) + +function rate_constant(V_app, m::Marcus, ::Val{true}; T = 298) + exp_arg = -(m.λ .- V_app).^2 ./ (4 .* m.λ .* kB * T) + m.A .* exp.(exp_arg) +end + +function rate_constant(V_app, m::Marcus, ::Val{false}; T = 298) + exp_arg = -(m.λ .+ V_app).^2 ./ (4 .* m.λ .* kB * T) + m.A .* exp.(exp_arg) +end diff --git a/src/kinetic_models/MarcusHushChidsey.jl b/src/kinetic_models/MarcusHushChidsey.jl new file mode 100644 index 0000000..6afd8f5 --- /dev/null +++ b/src/kinetic_models/MarcusHushChidsey.jl @@ -0,0 +1,36 @@ +""" + MarcusHushChidsey(A, λ, average_dos) + MarcusHushChidsey(λ, average_dos) + MarcusHushChidsey(λ) + +Computes Marcus-Hush-Chidsey kinetics: 10.1126/science.251.4996.919 + +Note that for "typical" reorganization energy values (in the vicinity of 10*kT at typical temperatures, i.e. a few tenths of an eV), `AsymptoticMarcusHushChidsey` is comparably accurate to and much faster to evaluate than this model. + +If either the prefactor or the average dos are omitted, their values are assumed to be 1. Note that strictly speaking, `average_dos` and the prefactor `A` are redundant. They are both included primarily to facilitate comparisons with similarly parametrized Marcus-like models such as `MarcusHushChidseyDOS`. +""" +struct MarcusHushChidsey{T} <: IntegralModel{T} + A::T + λ::T + average_dos::T + function MarcusHushChidsey(A, λ, average_dos) + ps = consistent_params(Float32.(A), Float32.(λ), Float32.(average_dos)) + new{typeof(ps[1])}(ps...) + end +end + +# default prefactor is 1 +MarcusHushChidsey(A, λ) = MarcusHushChidsey(A, λ, 1.0) +# assume prefactor = 1 and avg_dos = 1 +MarcusHushChidsey(λ) = MarcusHushChidsey(1.0, λ, 1.0) +# convert more detailed DOS information to just pull out average +MarcusHushChidsey(A, λ, dd::DOSData) = MarcusHushChidsey(A, λ, dd.average_value) +MarcusHushChidsey(A, λ, dos_file::Union{Matrix,String}; kwargs...) = + MarcusHushChidsey(A, λ, DOSData(dos_file; kwargs...)) + +# TODO: Check that both this and +DOS versions still match original paper +function integrand(mhc::MarcusHushChidsey, V_dl, ox::Val; T = 298) + marcus_term(E, ::Val{true}) = exp.(-((E .- mhc.λ .+ V_dl) .^ 2) ./ (4 * mhc.λ * kB * T)) + marcus_term(E, ::Val{false}) = exp.(-((E .- mhc.λ .- V_dl) .^2) ./ (4 * mhc.λ * kB * T)) + E -> mhc.A .* mhc.average_dos .* marcus_term(E, ox) .* fermi_dirac.(E; T=T) +end \ No newline at end of file diff --git a/src/kinetic_models/MarcusHushChidseyDOS.jl b/src/kinetic_models/MarcusHushChidseyDOS.jl new file mode 100644 index 0000000..cb8fca8 --- /dev/null +++ b/src/kinetic_models/MarcusHushChidseyDOS.jl @@ -0,0 +1,127 @@ +include("../utils/quantum_capacitance.jl") + +""" + MarcusHushChidseyDOS(A=1.0, λ, dos) + MarcusHushChidseyDOS(A=1.0, λ, dos_file) + +Computes Marcus-Hush-Chidsey + DOS kinetics as described in Kurchin and Viswanathan: 10.1063/5.0023611 + +NB: At the moment, this will allow for vector `A` and `λ` parameters, but will presume that all models correspond to the same DOS. +""" +struct MarcusHushChidseyDOS{T} <: IntegralModel{T} + A::T + λ::T + dos::DOSData + function MarcusHushChidseyDOS(A, λ, dos) + ps = consistent_params(Float32.(A), Float32.(λ)) + new{typeof(ps[1])}(ps..., dos) + end +end + +function Base.show(io::IO, mhcd::MarcusHushChidseyDOS) + s = repr(typeof(mhcd)) * "(" + s *= "A=$(round.(mhcd.A, sigdigits=3)), λ=$(round.(mhcd.λ, sigdigits=3)))" + print(io, s) +end + +# default prefactor is 1 +MarcusHushChidseyDOS(λ, dd::DOSData) = MarcusHushChidseyDOS(1.0, λ, dd) + +MarcusHushChidseyDOS(A, λ, dos_file::Union{Matrix,String}; kwargs...) = + MarcusHushChidseyDOS(A, λ, DOSData(dos_file; kwargs...)) + +MarcusHushChidseyDOS(λ, dos_file::Union{Matrix,String}; kwargs...) = MarcusHushChidseyDOS(1.0, λ, DOSData(dos_file; kwargs...)) + +# TODO: make one version of this for the whole package, currently there are some sign convention differences between this and the MHC version +marcus_term(V, E, λ, ::Val{true}, T) = exp.(-(( (λ .- V) .+ E) .^ 2) ./ (4 * λ * kB * T)) +marcus_term(V, E, λ, ::Val{false}, T) = exp.(-(( (λ .+ V) .- E) .^ 2) ./ (4 * λ * kB * T)) + +# Fermi-Dirac function or the complement thereof +fd(E, ::Val{true}, T) = 1 .- fermi_dirac(E; T = T) +fd(E, ::Val{false}, T) = fermi_dirac(E; T = T) + +function integrand( + mhcd::MarcusHushChidseyDOS, + V_dl, + ox::Val; + T = 298, + V_q = 0.0, +) + E -> mhcd.A .* mhcd.dos.interp_func.(E .+ V_q) .* fd(E, ox, T) .* marcus_term(V_dl, E, mhcd.λ, ox, T) +end + +function rate_constant( + V_app, + model::MarcusHushChidseyDOS, + args...; + T = 298, + E_min = model.dos.E_min, + E_max = model.dos.E_max, + calc_cq::Bool = false, + C_dl = 10.0, + Vq_min = -0.5, + Vq_max = 0.5, + kwargs..., +) + if calc_cq + rate_constant_cq( + V_app, + model, + args...; + C_dl = C_dl, + Vq_min = Vq_min, + Vq_max = Vq_max, + T = T, + E_min = min(E_min, E_min-Vq_min), + E_max = max(E_max, E_max-Vq_max), + ) + else + n, w = scale(E_min, E_max) + f = integrand(model, V_app, args...; T = T) + sum(w .* f.(n)) + end +end + +""" + rate_constant_cq(E_min, E_max, V_app, model::MarcusHushChidseyDOS, ox::Bool; C_dl=10.0, Vq_min=-0.5, Vq_max=0.5, kwargs...) + rate_constant_cq(E_min, E_max, V_app, model::MarcusHushChidseyDOS; C_dl=10.0, Vq_min=-0.5, Vq_max=0.5, kwargs...) + +Compute the rate constant k predicted by a `MarcusHushChidseyDOS` model at a applied voltage `V_app`, including the effects of quantum capacitance. 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. +""" +function rate_constant_cq( + V_app, + model::MarcusHushChidseyDOS, + ox::Bool; + C_dl = 10.0, + Vq_min = -0.5, + Vq_max = 0.5, + T = 298, + E_min = model.dos.E_min, + E_max = model.dos.E_max, +) + V_dl_interp = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl) + V_dl = V_dl_interp(V_app) + V_q = V_app - V_dl + n, w = scale(E_min, E_max) + f = integrand(model, V_dl, ox; T = T, V_q = V_q) + sum(w .* f.(n)) +end + +# TODO: clean these up...but add tests first +function rate_constant_cq( + V_app, + model::MarcusHushChidseyDOS; + C_dl = 10.0, + Vq_min = -0.5, + Vq_max = 0.5, + T = 298, + E_min = model.dos.E_min, + E_max = model.dos.E_max, +) + V_dl_interp = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl) + V_dl = V_dl_interp(V_app) + V_q = V_app - V_dl + n, w = scale(E_min, E_max) + f = integrand(model, V_dl; T = T, V_q = V_q) + sum(w .* f.(n)) +end diff --git a/src/kinetic_models/kinetic_models.jl b/src/kinetic_models/kinetic_models.jl new file mode 100644 index 0000000..2050cd3 --- /dev/null +++ b/src/kinetic_models/kinetic_models.jl @@ -0,0 +1,111 @@ +using .DOS +using SpecialFunctions +using Statistics +using Interpolations +using DelimitedFiles + +include("../utils/misc.jl") + +""" + fermi_dirac(E, T=298) + +Compute the value of the Fermi-Dirac distribution at energy `E` (relative to the Fermi energy) and temperature `T`. +""" +fermi_dirac(E; T = 298) = inv.(1 .+ exp.(E ./ (kB * T))) + +""" + KineticModel + +Abstract base class for kinetic models. +""" +abstract type KineticModel{T} end + +# return a new one with a scaled prefactor +import Base.* +function *(c, km::KineticModel) + new_A = c*km.A + other_params = getfield.([km], propertynames(km))[2:end] + eval(nameof(typeof(km)))(new_A, other_params...) +end + +import Base.length +length(km::KineticModel) = length(km.A) + +# generic pretty printing +function Base.show(io::IO, m::KineticModel) + s = repr(typeof(m)) * "(" + for field in propertynames(m) + s = + s * + string(field) * + "=" * + string(round.(getproperty(m, field), sigdigits = 3)) * + ", " + end + s = s[1:end-2] * ")" + print(io, s) +end + +""" + NonIntegralModel + +Abstract base class for kinetic models whose rates can be computed directly from an input voltage without requiring an energy integral. All subtypes must dispatch the `rate_constant` function. +""" +abstract type NonIntegralModel{T} <: KineticModel{T} end + +""" + IntegralModel + +Abstract base class for "Marcus-like" kinetic models that require computation of an energy integral. All subtypes need to dispatch the `rate_constant` function directly, or dispatch the `integrand` function and use the default `rate_constant` dispatch. +""" +abstract type IntegralModel{T} <: KineticModel{T} end + +# check to catch missed dispatches for new types +# 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 +# dispatch for net rates +integrand(km::IntegralModel, V; kwargs...) = + E -> integrand(km, V, Val(true); kwargs...)(E) - integrand(km, V, Val(false); kwargs...)(E) +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(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). + +If the model is an `IntegralModel`, integration bounds `E_min` and `E_max` may be supplied as kwargs. Integration is done via GK quadrature. + +If calc_cq flag is passed, optionally compute voltage shifts due to quantum capacitance (only applicable to `MarcusHushChidseyDOS` models). +""" +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...) + +# TODO: add tests that both args and kwargs are correctly captured here (also for the Val thing) +# "callable" syntax +(m::KineticModel)(V_app, args...; kwargs...) = rate_constant(V_app, m, args...; kwargs...) + +function rate_constant( + V_app, + model::IntegralModel, + args...; # would just be the ox flag, if present + T = 298, + E_min = -100 * kB * T, + E_max = 100 * kB * T +) + n, w = scale(E_min, E_max) + f = integrand(model, V_app, args...; T = T) + sum(w .* f.(n)) +end + +include("ButlerVolmer.jl") +include("Marcus.jl") +include("AsymptoticMarcusHushChidsey.jl") +include("MarcusHushChidsey.jl") +include("MarcusHushChidseyDOS.jl") diff --git a/src/phase_diagrams.jl b/src/phase_diagrams.jl index 0b6bddb..5acbe60 100644 --- a/src/phase_diagrams.jl +++ b/src/phase_diagrams.jl @@ -4,39 +4,38 @@ using NLsolve using FastGaussQuadrature # define some constants and default parameters -const kB = 8.617e-5 const room_T = 298 const muoA_default = 0.02 const muoB_default = 0.03 const Ω_default = 0.1 -N = 1000 -quadfun = gausslegendre # our familiar thermodynamic functions h(x;Ω=Ω_default) = @. x*(1-x)*Ω # enthalpy of mixing -s(x) = @. -kB*(x*log(x+eps(Float64)) + (1-x)*log(1-x+eps(Float64))) # entropy per particle...added epsilons in the hopes that things will be less obnoxious at the edges -g_thermo(x; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) = @. h(x;Ω=Ω_default) - T*s(x)+ muoA*(1-x) + muoB*x # Gibbs free energy per particle +s(x) = @. -kB*(x*log(x+eps(Float32)) + (1-x)*log(1-x+eps(Float32))) # entropy per particle...added epsilons in the hopes that things will be less obnoxious at the edges +g_thermo(x; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) = @. h(x;Ω) - T*s(x)+ muoA*(1-x) + muoB*x # Gibbs free energy per particle μ_thermo(x; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) = @. (1-2*x)*Ω + kB*T*log(x/(1-x)) + muoB-muoA # chemical potential +prefactor(x, intercalate::Bool) = intercalate ? (1 .- x) : x + """ µ 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). """ -function µ_kinetic(I, km::KineticModel; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) - thermo_term(x) = μ_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) - μ(x::Real) = thermo_term(x) .+ fit_overpotential((1-x)*km, I) - μ(x::AbstractVector) = thermo_term(x) .+ fit_overpotential((1 .- x).*Ref(km), I) +function µ_kinetic(I, km::KineticModel; intercalate=true, kwargs...) + thermo_term(x) = μ_thermo(x; kwargs...) + μ(x::Real) = thermo_term(x) .+ overpotential(I, prefactor(x, intercalate)*km) + μ(x::AbstractVector) = thermo_term(x) .+ overpotential(I, prefactor(x, intercalate).*Ref(km)) return μ end -function g_kinetic(I, km::KineticModel; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) - thermo_term(x) = g_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) - #TODO: gradient of this term is just value of fit_overpotential(x) +function g_kinetic(I, km::KineticModel; intercalate=true, kwargs...) + thermo_term(x) = g_thermo(x; kwargs...) + #TODO: gradient of this term is just value of overpotential(x) function kinetic_term(x) - f(x) = ElectrochemicalKinetics.fit_overpotential( (1 .- x) .* Ref(km), I) - n, w = ElectrochemicalKinetics.scale(zero.(x), x) + f(x) = ElectrochemicalKinetics.overpotential(I, prefactor(x, intercalate) * km) + n, w = ElectrochemicalKinetics.scale_coarse(zero.(x), x) map((w, n) -> sum(w .* f(n)), eachcol(w), eachcol(n)) end g(x) = thermo_term(x) .+ kinetic_term(vec(x)) @@ -46,46 +45,45 @@ end # 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; Ω=Ω_default, muoA=muoA_default, muoB=muoB_default, T=room_T) - g = g_kinetic(I, km; Ω=Ω, muoA=muoA, muoB=muoB, T=T) - µ = µ_kinetic(I, km; Ω=Ω, muoA=muoA, muoB=muoB, T=T) +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...) [(g(x[2]) - g(x[1]))/(x[2] - x[1]) .- μ(x[1]), μ(x[2])-μ(x[1])] end -# case where we want to check many points at once (shape of x should be N x 2) -function common_tangent(x::Array, I, km::KineticModel; kwargs...) - g = g_kinetic(I, km; kwargs...) - µ = µ_kinetic(I, km; kwargs...) - Δg = g(x[:,2]) - g(x[:,1]) - Δx = x[:,2] - x[:,1] - μ1 = μ(x[:,1]) - Δμ = μ(x[:,2]) - μ1 - hcat(Δg ./ Δx .- μ1, Δμ) -end - -function find_phase_boundaries(I, km::KineticModel; quadfun=quadfun, N=N, kwargs...) +# 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 myct!(storage, x) - res = common_tangent(x, I, km; kwargs...) + res = common_tangent(x, I, km; intercalate=intercalate, kwargs...) storage[1] = res[1] storage[2] = res[2] end - # TODO: grad! here from AD - x1 = nlsolve(myct!, [0.05 0.95]) + x1 = nlsolve(myct!, guess, show_trace=verbose, ftol=1e-6, xtol=1e-6) x1.zero end -function phase_diagram(km::KineticModel; kwargs...) - # TODO: write this, lol +# TODO: add T_max to handle low-T BV cases +function phase_diagram(km::KineticModel; I_step=10, verbose=false, intercalate=true, kwargs...) + I = 0 + pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, kwargs...) + pbs = pbs_here' + I_vals = [0] + while abs(pbs_here[2] - pbs_here[1]) > 1e-3 + I = I + I_step + if verbose + println("Solving at I=", I, "...") + end + try + pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_here, kwargs...) + pbs = vcat(pbs, pbs_here') + push!(I_vals, I) + if verbose + println("Phase boundaries are at ", pbs_here) + end + catch e + println("Solve failed at I=", I) + end + end + return vcat(pbs[:,1], reverse(pbs[:,2])), vcat(I_vals, reverse(I_vals)) end - -# bv = ButlerVolmer(300) -# I_vals = 10 .^ (1.1:0.025:3.1) - - -# this line takes a few seconds with B-V but aaages with anything else... -# pbs = [find_phase_boundaries(I, bv, T=330) for I in I_vals] - -# plot(vcat(pbs...), hcat(I_vals, I_vals), label="phase boundary") -# xlabel!("x") -# ylabel!("I") diff --git a/src/plot_fcns.jl b/src/plot_fcns.jl index 1381f3e..60305cf 100644 --- a/src/plot_fcns.jl +++ b/src/plot_fcns.jl @@ -9,7 +9,7 @@ Plot predicted rate constants for each model in the provided list. # Keyword Arguments * V_min and V_max: bounds of voltage, defaults to +/- 1.0 * plot_title -* kwargs for compute_k function +* kwargs for rate_constant function """ function plot_models( models::Vector{<:KineticModel}; @@ -20,7 +20,7 @@ function plot_models( ) V_range = range(V_min, V_max, length = 200) xs = repeat([V_range], length(models)) - ys = map(model -> [compute_k(V, model; kwargs...) for V in V_range], models) + ys = map(model -> abs.(rate_constant(V_range, model; kwargs...)), models) plot( xs, ys, @@ -42,7 +42,7 @@ Plot predicted rate constants for each model in the provided list. # Keyword Arguments * V_min and V_max: bounds of voltage, defaults to +/- 1.0 * plot_title -* kwargs for compute_k function +* kwargs for rate_constant function """ function plot_exp_and_models( exp_data::Matrix, @@ -56,7 +56,7 @@ function plot_exp_and_models( xs = Vector[V, repeat([V_range], length(models))...] ys = Vector[ exp_data[:, 2], - map(model -> compute_k(V_range, model; kwargs...), models)..., + map(model -> abs.(rate_constant(V_range, model; kwargs...)), models)..., ] # scatter plot of experimental data, lines for fits diff --git a/src/rate_constant.jl b/src/rate_constant.jl deleted file mode 100644 index 9ceef35..0000000 --- a/src/rate_constant.jl +++ /dev/null @@ -1,148 +0,0 @@ -using Statistics -using Interpolations -using DelimitedFiles - -# TODO: think about reaction direction and quantum cap dispatches; there's probably a cleaner way to handle them... - -""" - compute_k(V_app, model::KineticModel, ox::Bool; kwargs...) - compute_k(V_app, model::KineticModel; kwargs...) - compute_k(E_min, E_max, V_app, model::IntegralModel, ox::Bool; kwargs...) - compute_k(E_min, E_max, V_app, model::IntegralModel; kwargs...) - compute_k(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). - -If the model is an `IntegralModel`, integration bounds `E_min` and `E_max` must be supplied. Integration is done via GK quadrature. - -If calc_cq flag is passed, optionally compute voltage shifts due to quantum capacitance. -""" -compute_k(V_app, model::KineticModel, args...; kT = 0.026) = model(V_app, args...; kT = kT) - -compute_k( - V_app, - model::IntegralModel, - args...; # would just be the ox flag - kT = 0.026, - E_min = -100 * kT, - E_max = 100 * kT, - kwargs..., -) = begin - n, w = scale(E_min, E_max) - f = integrand(model, V_app, args...; kT = kT) - sum(w .* f.(n)) -end - -function compute_k( - V_app, - model::MarcusHushChidseyDOS, - ox::Bool; - kT = 0.026, - E_min = model.dos.E_min, - E_max = model.dos.E_max, - calc_cq::Bool = false, - C_dl = 10.0, - Vq_min = -0.5, - Vq_max = 0.5, - kwargs..., -) - if calc_cq - compute_k_cq( - V_app, - model, - ox; - C_dl = C_dl, - Vq_min = Vq_min, - Vq_max = Vq_max, - kT = kT, - E_min = min(E_min, E_min-Vq_min), - E_max = max(E_max, E_max-Vq_max), - ) - else - n, w = scale(E_min, E_max) - f = integrand(model, V_app, ox; kT = kT) - sum(w .* f.(n)) - end -end - -function compute_k( - V_app, - model::MarcusHushChidseyDOS; - kT = 0.026, - E_min = model.dos.E_min, - E_max = model.dos.E_max, - calc_cq::Bool = false, - C_dl = 10.0, - Vq_min = -0.5, - Vq_max = 0.5, - kwargs..., -) - if calc_cq - compute_k_cq( - V_app, - model; - kT = kT, - E_min = max(E_min, E_min-Vq_min), - E_max = min(E_max, E_max-Vq_max), - C_dl = C_dl, - Vq_min = Vq_min, - Vq_max = Vq_max, - ) - else - # invoke( - # compute_k, - # Tuple{typeof(V_app),IntegralModel}, - # V_app, - # model; - # kT = kT, - # E_min = E_min, - # E_max = E_max, - # ) - n, w = scale(E_min, E_max) - f = integrand(model, V_app; kT = kT) - sum(w .* f.(n)) - end -end - -""" - compute_k_cq(E_min, E_max, V_app, model::MarcusHushChidseyDOS, ox::Bool; C_dl=10.0, Vq_min=-0.5, Vq_max=0.5, kwargs...) - compute_k_cq(E_min, E_max, V_app, model::MarcusHushChidseyDOS; C_dl=10.0, Vq_min=-0.5, Vq_max=0.5, kwargs...) - -Compute the rate constant k predicted by a `MarcusHushChidseyDOS` model at a applied voltage `V_app`, including the effects of quantum capacitance. 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. -""" -function compute_k_cq( - V_app, - model::MarcusHushChidseyDOS, - ox::Bool; - C_dl = 10.0, - Vq_min = -0.5, - Vq_max = 0.5, - kT = 0.026, - E_min = model.dos.E_min, - E_max = model.dos.E_max, -) - V_dl_interp = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl) - V_dl = V_dl_interp(V_app) - V_q = V_app - V_dl - n, w = scale(E_min, E_max) - f = integrand(model, V_dl, ox; kT = kT, V_q = V_q) - sum(w .* f.(n)) -end - -function compute_k_cq( - V_app, - model::MarcusHushChidseyDOS; - C_dl = 10.0, - Vq_min = -0.5, - Vq_max = 0.5, - kT = 0.026, - E_min = model.dos.E_min, - E_max = model.dos.E_max, -) - V_dl_interp = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl) - V_dl = V_dl_interp(V_app) - V_q = V_app - V_dl - n, w = scale(E_min, E_max) - f = integrand(model, V_dl; kT = kT, V_q = V_q) - sum(w .* f.(n)) -end diff --git a/src/dos.jl b/src/utils/dos.jl similarity index 91% rename from src/dos.jl rename to src/utils/dos.jl index b947375..d5f7e9d 100644 --- a/src/dos.jl +++ b/src/utils/dos.jl @@ -3,7 +3,6 @@ module DOS using Interpolations using DelimitedFiles using Statistics -# using Zygote export DOSData, get_dos @@ -19,9 +18,9 @@ Struct for storing density of states (DOS) information. Can be constructed from """ struct DOSData interp_func - average_value::Float64 - E_min::Float64 - E_max::Float64 + average_value::Float32 + E_min::Float32 + E_max::Float32 end DOSData(dos_file; Ef=0, cut_energy=false) = DOSData(get_dos(dos_file; Ef, cut_energy)...) @@ -31,6 +30,9 @@ Base.show(io::IO, dd::DOSData) = print(io, "DOSData: avg value $(round(dd.averag (dd::DOSData)(E::Real) = dd.interp_func(E) +import Base.length +length(::DOSData) = 1 + """ get_dos(dos_file; Ef=0, cut_energy=false) @@ -44,10 +46,10 @@ Returns a callable interpolated DOS, the average value of DOS, and the lower and """ function get_dos(dos_file::String; kwargs...) # dos_data = Zygote.ignore() do - # x = readdlm(dos_file, Float64, skipstart=1) + # x = readdlm(dos_file, Float32, skipstart=1) # x # end - dos_data = readdlm(dos_file, Float64, skipstart=1) + dos_data = readdlm(dos_file, Float32, skipstart=1) get_dos(dos_data; kwargs...) end function get_dos(dd::Matrix; Ef=0, cut_energy=false) diff --git a/src/integration_utils.jl b/src/utils/integration_utils.jl similarity index 61% rename from src/integration_utils.jl rename to src/utils/integration_utils.jl index 43995bb..925910d 100644 --- a/src/integration_utils.jl +++ b/src/utils/integration_utils.jl @@ -16,14 +16,5 @@ function setup_integration(f, N) end end -const scale = setup_integration(gausslegendre, 1000) - -# Can't get this to broadcast correctly - but this is the right sort of -# approach to take here. Of note - this is 2x faster in the forwards pass -# and about 10x faster in the backwards pass compared to the other method -# function scale_integration_nodes2(unscaled_x, unscaled_w, lb, ub) -# w_scale = @. 0.5 * (ub - lb) -# x_scale = @. 0.5 * (ub + lb) + 0.5 * (ub - lb) -# unscaled_x * x_scale', unscaled_w * w_scale' -# # x_scale * unscaled_x, unscaled_w * w_scale' -# end +const scale = setup_integration(gausslegendre, 500) +const scale_coarse = setup_integration(gausslegendre, 50) diff --git a/src/utils/misc.jl b/src/utils/misc.jl new file mode 100644 index 0000000..370f106 --- /dev/null +++ b/src/utils/misc.jl @@ -0,0 +1,20 @@ +# make a list of parameters into a list of vectors of the same length +function consistent_params(ps...) + lengths = Set(length.(ps)) + if length(lengths) == 1 # they're already the same length + return ps + elseif (length(lengths)==2 && 1 in lengths) + len = [l for l in lengths if l!=1][1] + new_ps = [] + for p in ps + if length(p) == len + push!(new_ps, p) + else + push!(new_ps, repeat([p], len)) + end + end + return new_ps + else + throw(DimensionMismatch("Mixed length fields...they must be all the same length, or a combination of one length and length one ;)")) + end +end diff --git a/src/quantum_capacitance.jl b/src/utils/quantum_capacitance.jl similarity index 100% rename from src/quantum_capacitance.jl rename to src/utils/quantum_capacitance.jl diff --git a/test/fitting_tests.jl b/test/fitting_tests.jl index c2bd3dc..4b43c85 100644 --- a/test/fitting_tests.jl +++ b/test/fitting_tests.jl @@ -1,3 +1,4 @@ +# TODO: add tests for shortcutting in overpotential solve @testset "Overpotential Fitting" begin params = Dict(ButlerVolmer=>[], Marcus=>[0.25], AsymptoticMarcusHushChidsey=>[0.25], MarcusHushChidsey=>[0.25]) for model_type in [ButlerVolmer, Marcus, AsymptoticMarcusHushChidsey, MarcusHushChidsey] @@ -5,20 +6,20 @@ m = model_type(params[model_type]...) # test a number less than the initial guess V1 = 0.05 - k1 = compute_k(V1, m) - @test isapprox(fit_overpotential(m, k1)[1], V1, atol=1e-5) + k1 = rate_constant(V1, m) + @test isapprox(overpotential(k1, m), V1, atol=1e-5) # and one greater V2 = 0.15 - k2 = compute_k(V2, m) - @test isapprox(fit_overpotential(m, k2)[1], V2, atol=1e-5) - # test that vectorized fitting works...first for one model and multiple k's - # this test is turned off for now because better adjoint of fit_overpotential currently doesn't work with it... - # @test isapprox(fit_overpotential(m, [k1, k2]), [V1, V2], atol=5e-4) - # TODO: figure out sensible and consistent way to vectorize over k and/or V (e.g. use a matrix?) - # now test for one k and multiple models - @test isapprox(fit_overpotential([m, m], k1), [V1, V1], atol=5e-4) + k2 = rate_constant(V2, m) + @test isapprox(overpotential(k2, m), V2, atol=1e-5) + # test that vectorized fitting works, both for scalar k and vector model and vice versa + # TODO: figure out why this doesn't work for integral models + if model_type != MarcusHushChidsey + @test isapprox(overpotential([k1, k2], m), [V1, V2], atol=5e-4) + @test isapprox(overpotential(k1, [1,1]*m), [V1, V1], atol=5e-4) + end # make sure that both fails since indexing would be ambiguous - @test_throws MethodError fit_overpotential([m, m], [k1, k2]) + @test_throws MethodError overpotential([k1, k2], [m, m]) end end # for now, writing these separately to test different input DOSes @@ -33,13 +34,13 @@ for dos in [flat_dos, smooth_dos] m = MarcusHushChidseyDOS(0.25, dos) V1 = 0.05 - k1 = compute_k(V1, m) - @test isapprox(fit_overpotential(m, k1)[1], V1, atol=1e-5) + k1 = rate_constant(V1, m) + @test isapprox(overpotential(k1, m)[1], V1, atol=1e-5) V2 = 0.15 - k2 = compute_k(V2, m) - @test isapprox(fit_overpotential(m, k2)[1], V2, atol=1e-5) + k2 = rate_constant(V2, m) + @test isapprox(overpotential(k2, m)[1], V2, atol=1e-5) # this one is disabled for now for same reason as above - # @test isapprox(fit_overpotential(m, [k1, k2]), [V1, V2], atol=5e-4) + # @test isapprox(overpotential([k1, k2], m), [V1, V2], atol=5e-4) end end end diff --git a/test/integral_model_tests.jl b/test/integral_model_tests.jl index cab3ead..87a1ecb 100644 --- a/test/integral_model_tests.jl +++ b/test/integral_model_tests.jl @@ -1,31 +1,44 @@ # all these numbers are just references evaluated as of 2/22/22 - +# TODO: add tests of `integrand` function +# TODO: throughout, add tests of varying kT @testset "Integral Models" begin - @testset "MHC" begin + @testset "MarcusHushChidsey" begin mhc = MarcusHushChidsey(0.25) + @testset "Scalars" begin + # test some values + @test rate_constant(0, mhc) == 0.0 + @test rate_constant(0, mhc, true) == rate_constant(0, mhc, false) ≈ 0.005892741f0 + @test isapprox(rate_constant(0.1, mhc), 0.03062463, atol=1e-6) + + # test net rates and symmetry + @test isapprox(rate_constant(0.1, mhc, true) - rate_constant(0.1, mhc, false), rate_constant(0.1, mhc), atol=1e-6) + @test -rate_constant(-0.1, mhc) == rate_constant(0.1, mhc) + end - # test some values - @test compute_k(0, mhc) == 0.0 - @test compute_k(0, mhc, true) == compute_k(0, mhc, false) ≈ 0.0061371322 - @test isapprox(compute_k(0.1, mhc), 0.031239978, atol=1e-6) - - # test net rates and symmetry - @test isapprox(compute_k(0.1, mhc, true) - compute_k(0.1, mhc, false), compute_k(0.1, mhc), atol=1e-6) - @test compute_k(-0.1, mhc) == compute_k(0.1, mhc) + @testset "Vector Voltages" begin + test_vector_voltages(mhc, Vs) + @test isapprox(mhc(Vs), -mhc(-Vs)) #symmetry + # test that it's close to asymptotic version + amhc = AsymptoticMarcusHushChidsey(0.25) + @test all(isapprox.(rate_constant(Vs, mhc), rate_constant(Vs, amhc), atol=2e-4)) + end - # test that it's close to asymptotic version, also that vector inputs work appropriately - amhc = AsymptoticMarcusHushChidsey(0.25) - @test all(isapprox.(compute_k(Vs, mhc), compute_k(Vs, amhc), atol=2e-4)) + @testset "Vector Models" begin + A_vals = [1, 10, 1000] + λ_vals = [0.1, 0.4, 0.7] + avg_dos_val = 1.0 + test_vector_models(MarcusHushChidsey, [A_vals, λ_vals, avg_dos_val]) + end end - @testset "MHC+DOS" begin + @testset "MarcusHushChidseyDOS" begin # test that a uniform DOS version matches MHC flat_dos = [-5 1; 5 1] mhcd = MarcusHushChidseyDOS(0.25, flat_dos) - @test isapprox(compute_k(0, mhcd), 0.0, atol=1e-6) + @test isapprox(rate_constant(0, mhcd), 0.0, atol=1e-6) mhc = MarcusHushChidsey(0.25) - mhcd_k = compute_k(Vs, mhcd) - mhc_k = compute_k(Vs, mhc) + mhcd_k = rate_constant(Vs, mhcd) + mhc_k = rate_constant(Vs, mhc) @test all(isapprox.(mhc_k, mhcd_k, atol=1e-6)) # make some cartoon smooth DOSes... @@ -35,11 +48,14 @@ dos = hcat(Es, gaussian.(0.7 .*(Es .- 2)) .+ gaussian.(0.7 .*(Es .+ 2))) mhcd = MarcusHushChidseyDOS(0.25, dos) # this should be symmetric - @test all(isapprox.(compute_k(Vs, mhcd), compute_k(-Vs, mhcd), atol=1e-6)) + @test all(isapprox.(rate_constant(Vs, mhcd), -rate_constant(-Vs, mhcd), atol=1e-6)) # now try an asymmetric one dos = hcat(Es .+ 1, dos[:,2]) mhcd = MarcusHushChidseyDOS(0.25, dos) - @test all(compute_k(-Vs, mhcd) .> compute_k(Vs, mhcd)) + @test all(-rate_constant(-Vs, mhcd) .> rate_constant(Vs, mhcd)) + + params = [[2, 20], [0.05, 0.1], DOSData(dos)] + test_vector_models(MarcusHushChidseyDOS, params) end @testset "Quantum Capacitance" begin diff --git a/test/non_integral_model_tests.jl b/test/non_integral_model_tests.jl index 0f970a1..49f3166 100644 --- a/test/non_integral_model_tests.jl +++ b/test/non_integral_model_tests.jl @@ -1,61 +1,94 @@ # test values come either from analytical inversion of expression (BV) or references evaluated as of 2/22/22 and presumed to be correct +# TODO: throughout, add tests of varying kT + @testset "Non-Integral Models" begin @testset "Butler-Volmer" begin bv1 = ButlerVolmer() - bv2 = ButlerVolmer(3.0) - bv3 = ButlerVolmer(1.0, 0.2) + bv2 = ButlerVolmer(3.0, 0.5) + bv3 = ButlerVolmer(0.2) - @testset "Scalar Arguments" begin + @testset "Scalars" begin # zero voltage - @test all(compute_k.(Ref(0), [bv1, bv2, bv3]).==0) + @test all(rate_constant.(Ref(0), [bv1, bv2, bv3]).==0) # oxidative for bv in [bv1, bv2, bv3] - @test bv(.026/bv.α, Val(true)) ≈ bv.A * ℯ + @test bv(kB*298/bv.α, true) ≈ bv.A * ℯ end # reductive for bv in [bv1, bv2, bv3] - @test bv(-.026/(1-bv.α), Val(false)) ≈ bv.A * ℯ + @test bv(-kB*298/(1-bv.α), false) ≈ bv.A * ℯ end # net for bv in [bv1, bv2, bv3] - @test bv(.026/bv.α) ≈ bv.A * (ℯ - exp((bv.α-1)/bv.α)) + @test bv(kB*298/bv.α) ≈ bv.A * (ℯ - exp((bv.α-1)/bv.α)) end end - @testset "Vector Arguments" begin - @test bv1(.026 .* collect(0:2:6), Val(true)) == exp.(0:3) - @test bv1(.026 .* collect(0:2:6)) == exp.(0:3) - exp.(-(0:3)) - # TODO: test bv2, bv3 + @testset "Vector Voltages" begin + test_vector_voltages(bv1, Vs) + test_vector_voltages(bv2, Vs) + test_vector_voltages(bv3, Vs) + for model in [bv1, bv2] # bv3 won't be symmetrical + @test isapprox(model(Vs), -model(-Vs)) #symmetry + end + end + + @testset "Vector Models" begin + A_vals = [0.7, 2.0] + α_vals = [0.4, 0.5, 0.6] + @test_throws DimensionMismatch ButlerVolmer(A_vals, α_vals) # invalid, different lengths + + test_vector_models(ButlerVolmer, [A_vals, α_vals[1:2]]) end end @testset "Marcus" begin - # TODO: this could all be more thorough, I literally just picked some arbitrary points - m = Marcus(0.25) - @test m(0.0) == 0.0 - @test isapprox(m(0.05, Val(true)), 0.031381, atol=1e-6) - @test isapprox(m(0.2, Val(false)), 0.908324, atol=1e-6) - @test isapprox(m(0.13), 0.570862, atol=1e-6) - @test m(0.13) == m(-0.13) - - # test vector voltages and also symmetry - @test m(Vs) == m(-Vs) + m1 = Marcus(0.25) + m2 = Marcus(10, 0.25) + @testset "Scalars" begin + # TODO: this could all be more thorough, I literally just picked some arbitrary points + @test m1(0.0) == 0.0 + @test isapprox(m1(0.05, false), 0.030055, atol=1e-6) + @test isapprox(m1(0.2, true), 0.907235, atol=1e-6) + @test isapprox(m1(0.13), 0.5671645, atol=1e-6) + @test isapprox(m2(0.13), 5.671645, atol=1e-6) + end + + @testset "Vector Voltages" begin + test_vector_voltages(m1, Vs) + test_vector_voltages(m2, Vs) + for model in [m1, m2] + @test isapprox(model(Vs), -model(-Vs)) #symmetry + end + end + @testset "Vector Models" begin + A_vals = [1, 100, 10000] + λ_vals = [0.1, 0.3, 0.7] + test_vector_models(Marcus, [A_vals, λ_vals]) + end end @testset "asymptotic MHC" begin amhc = AsymptoticMarcusHushChidsey(0.25) - @test amhc(0.0) == 0.0 - @test amhc(0.0, true) == amhc(0.0, false) ≈ 0.00596440064 - @test amhc(0.2) == amhc(-0.2) ≈ 0.1006330238 - @test amhc(0.2, kT=.015) ≈ 0.06360960459 - @test amhc(0.2, kT=.03) ≈ 0.11255129857 - - vec_out = amhc(Vs) - @test length(vec_out) == length(Vs) - @test vec_out[1] == amhc(Vs[1]) + @testset "Scalars" begin + @test amhc(0.0) == 0.0 + @test amhc(0.0, true) == amhc(0.0, false) ≈ 0.00573491969678f0 + @test amhc(0.2) == -amhc(-0.2) ≈ 0.09964746f0 + @test amhc(0.2, T=175) ≈ 0.063909331f0 + @test amhc(0.2, T=350) ≈ 0.11301532f0 + end + @testset "Vector Voltages" begin + test_vector_voltages(amhc, Vs) + @test isapprox(amhc(Vs), -amhc(-Vs)) #symmetry + end + @testset "Vector Models" begin + A_vals = [1, 100, 10000] + λ_vals = [0.1, 0.3, 0.7] + test_vector_models(AsymptoticMarcusHushChidsey, [A_vals, λ_vals]) + end end end diff --git a/test/phase_diagram_tests.jl b/test/phase_diagram_tests.jl index 2f4ce17..37a6e70 100644 --- a/test/phase_diagram_tests.jl +++ b/test/phase_diagram_tests.jl @@ -4,35 +4,38 @@ muoA = 0.02 muoB = 0.03 T = 298 +# TODO: add tests where we change Ω, etc. +# TODO: add tests for deintercalation case + @testset "Basic free energy functions" begin # enthalpy of mixing @test h(0.5) == 0.025 # value @test h(0.1) ≈ h(0.9) # symmetry - @test h([0.1, 0.2, 0.7]) ≈ [0.009, 0.016, 0.021] # vector input + @test all(h([0.1, 0.2, 0.7]) .≈ [0.009, 0.016, 0.021]) # vector input @test h(0.5; Ω=1) == 0.25 # kwargs # similar set for entropy... - @test s(0.5) ≈ -kB * 2 * 0.5 * log(0.5) + @test isapprox(s(0.5f0), -kB * log(0.5f0), atol=1f-5) @test s(0.1) == s(0.9) - @test s([0.1, 0.2, 0.7]) ≈ [2.8012399817e-5, 4.3119676836e-5, 5.2638176908e-5] + @test all(s([0.1, 0.2, 0.7]) .≈ [2.8013f-5, 4.3121f-5, 5.264f-5]) # ...and thermodynamic free energy... - @test g_thermo(0.5) == h(0.5) - T * s(0.5) + muoA * 0.5 + muoB * 0.5 ≈ 0.032200909 - @test g_thermo(0.0) ≈ muoA - @test g_thermo(1.0) ≈ muoB - @test g_thermo([0.1, 0.2, 0.7]) ≈ [0.02165230485, 0.0251503363, 0.032313823] + @test g_thermo(0.5) == h(0.5) - T * s(0.5) + muoA * 0.5 + muoB * 0.5 ≈ 0.0322002f0 + @test isapprox(g_thermo(0.0), muoA, atol=1e-8) + @test isapprox(g_thermo(1.0), muoB, atol=1e-8) + @test all(g_thermo([0.1, 0.2, 0.7]) .≈ [0.021651988f0, 0.02514984585595f0, 0.0323132232f0]) # ...and thermodynamic chemical potential. @test all(isinf.(μ_thermo([0.0, 1.0]))) @test μ_thermo(0.5) == μ_thermo(0.5, T=500) ≈ muoB - muoA - @test μ_thermo(0.1) ≈ 0.033578217 - @test μ_thermo(0.9, T=400) ≈ 0.0057339367 - @test μ_thermo([0.1, 0.2, 0.7]) ≈ [0.033578217, 0.034401818, -0.008242526] - @test μ_thermo([0.1, 0.2, 0.7], T=350) ≈ [0.023732805, 0.028190055, -0.00444592] + @test μ_thermo(0.1) ≈ 0.033576035f0 + @test isapprox(μ_thermo(0.9, T=400), 0.00573686571, atol=5e-6) + @test all(μ_thermo([0.1, 0.2, 0.7]) .≈ [0.033576f0, 0.03440044f0, -0.0082417f0]) + @test all(μ_thermo([0.1, 0.2, 0.7], T=350) .≈ [0.02373024f0, 0.0281884f0, -0.00444493f0]) end # prefactors are set so that k vs. η plots roughly line up for small η -bv = ButlerVolmer(300) +bv = ButlerVolmer(300, 0.5) m = Marcus(5000, 0.3) amhc = AsymptoticMarcusHushChidsey(70000, 0.3) kms = [bv, m, amhc] @@ -42,28 +45,28 @@ xs = [0.1, 0.5, 0.95] @testset "Kinetic μ" begin # all these numbers are just references evaluated as of 2/24/22 μ_50_vals = Dict( - ButlerVolmer=>Dict( - 0.1 => 0.0383864736, - 0.5 => 0.01862765755, - 0.95 => 0.06237022288), - Marcus=>Dict( - 0.1 => 0.02841242588577, - 0.5 => 0.01928206795, - 0.95 => 0.07492901629449), - AsymptoticMarcusHushChidsey=>Dict( - 0.1 => 0.0389118763, - 0.5 => 0.0195713281175, - 0.95 => 0.0735097011) + :ButlerVolmer=>Dict( + 0.1 => 0.038325055f0, + 0.5 => 0.018521359f0, + 0.95 => 0.0615507f0), + :Marcus=>Dict( + 0.1 => 0.03886537f0, + 0.5 => 0.01950138f0, + 0.95 => 0.0760257f0), + :AsymptoticMarcusHushChidsey=>Dict( + 0.1 => 0.0390861f0, + 0.5 => 0.01988688f0, + 0.95 => 0.07516774f0) ) μ_200_vals = Dict( - ButlerVolmer=>[0.0524237924, 0.04250974, 0.13058880458], - Marcus=>[0.054009696387679, 0.045881168639779, 0.212205862750], - AsymptoticMarcusHushChidsey=>[0.05451269755, 0.046343172588, 0.17452291988558] + :ButlerVolmer=> [0.0521894f0, 0.0421092f0, 0.1289288f0], + :Marcus=> [0.0544722f0, 0.0466285f0, 0.2127188f0], + :AsymptoticMarcusHushChidsey=>[0.055173985f0, 0.04740624957f0, 0.17594977f0] ) μ_100_T400_vals = Dict( - ButlerVolmer=>[0.02384219096, 0.02702875, 0.1212749347], - Marcus=>[0.024573427974, 0.028429814254, 0.1529930329757], - AsymptoticMarcusHushChidsey=>[0.024890385314, 0.028908839244, 0.14468140815] + :ButlerVolmer=>[0.023721278f0, 0.02681895f0, 0.120048858f0] , + :Marcus=>[0.02481338f0, 0.0288534856f0, 0.1539812383f0], + :AsymptoticMarcusHushChidsey=>[0.0252367557f0, 0.029513253f0, 0.146351765f0] ) for km in kms @testset "$(typeof(km))" begin @@ -71,25 +74,25 @@ xs = [0.1, 0.5, 0.95] μ_200 = μ_kinetic(200, km) μ_100_T400 = μ_kinetic(100, km, T=400) # test that the right (approximate) relationships hold - μs = μ_50_vals[typeof(km)] + μs = μ_50_vals[typeof(km).name.name] for x in xs - @test μ_50(x) ≈ μs[x] ≈ μ_thermo(x) + fit_overpotential((1-x)*km, 50) + @test μ_50(x) ≈ μs[x] ≈ μ_thermo(x) + overpotential(50, (1-x)*km) end # test vector inputs - @test μ_200(xs) ≈ μ_200_vals[typeof(km)] - @test μ_100_T400(xs) == μ_100_T400.(xs) ≈ μ_100_T400_vals[typeof(km)] + @test μ_200(xs) ≈ μ_200_vals[typeof(km).name.name] + @test μ_100_T400(xs) == μ_100_T400.(xs) && all(isapprox.(µ_100_T400(xs), μ_100_T400_vals[typeof(km).name.name], atol=1f-5)) end end end @testset "Kinetic g" begin g_200_T400_vals = Dict( - ButlerVolmer => [0.0205857425, 0.037693617, 0.06661696], - Marcus => [0.02073470111, 0.03874647481, 0.07399607032], - AsymptoticMarcusHushChidsey => [0.0207838185, 0.0390031065, 0.0729910686] + :ButlerVolmer => [0.020563, 0.03755, 0.0661336], + :Marcus => [0.0207786, 0.0390248, 0.074701], + :AsymptoticMarcusHushChidsey => [0.0208467, 0.03939946, 0.0740394] ) - g_50_2_vals = Dict(ButlerVolmer=>0.03519728018258, Marcus=>0.03395267141265, AsymptoticMarcusHushChidsey=>0.0347968044) + g_50_2_vals = Dict(:ButlerVolmer=>0.03515968, :Marcus=>0.035497, :AsymptoticMarcusHushChidsey=>0.0356339) for km in kms @testset "$(typeof(km))" begin μ_50 = μ_kinetic(50, km) @@ -102,10 +105,10 @@ end # check scalar input works @test g_50(xs[2]) == g_50(xs)[2] - @test isapprox(g_50(xs)[2], g_50_2_vals[typeof(km)], atol=1e-5) + @test isapprox(g_50(xs)[2], g_50_2_vals[typeof(km).name.name], atol=1e-5) # check a few other values - @test all(isapprox.(g_200_T400(xs), g_200_T400_vals[typeof(km)], atol=1e-5)) + @test all(isapprox.(g_200_T400(xs), g_200_T400_vals[typeof(km).name.name], atol=1e-5)) end end end @@ -128,10 +131,7 @@ end @test v3[2] < v1[2] # test actual numerical values too - @test isapprox(v1, [0.04584 0.839942], atol=1e-6) - @test isapprox(v2, [0.0837787 0.795035], atol=1e-6) - @test isapprox(v3, [0.107585 0.675948], atol=1e-6) - - # TODO: next for multiple currents at once (may require some syntax tweaks) - + @test all(isapprox.(v1, [0.045547, 0.841501], atol=1e-5)) + @test all(isapprox.(v2, [0.083289, 0.796802], atol=1e-4)) + @test all(isapprox.(v3, [0.105747, 0.6809], atol=1e-4)) end diff --git a/test/runtests.jl b/test/runtests.jl index 0630324..9015c45 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,28 @@ using ElectrochemicalKinetics using QuadGK using Test +function test_vector_voltages(model, voltages) + @test isapprox(rate_constant(voltages, model), rate_constant.(voltages, Ref(model))) + @test isapprox(rate_constant(voltages, model, T=1000), rate_constant.(voltages, Ref(model), T=1000)) +end + +function test_vector_models(ModelType, params) + model1 = ModelType(params[2:end]...) + model2 = ModelType(params...) + broadcast_params = [] + for p in params + if length(p) == 1 + push!(broadcast_params, Ref(p)) + else + push!(broadcast_params, p) + end + end + @test isapprox(rate_constant(0.1, model1), rate_constant.(Ref(0.1), ModelType.(broadcast_params[2:end]...))) + @test isapprox(rate_constant(0.1, model1, T=1000), rate_constant.(Ref(0.1), ModelType.(broadcast_params[2:end]...), T=1000)) + @test isapprox(rate_constant(0.1, model2), rate_constant.(Ref(0.1), ModelType.(broadcast_params...))) + @test isapprox(rate_constant(0.1, model2, T=200), rate_constant.(Ref(0.1), ModelType.(broadcast_params...), T=200)) +end + @testset "ElectrochemicalKinetics.jl" begin global Vs = 0.01:0.01:0.1 @testset "Rate constant computation" begin