diff --git a/.gitignore b/.gitignore index b067edd..0acd06b 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ -/Manifest.toml +./Manifest.toml +./github diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 0000000..14d21c1 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,303 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.11.1" +manifest_format = "2.0" +project_hash = "0978e35ac3c4ed2f2567f101098e29d371b1f36c" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +version = "1.11.0" + +[[deps.Chain]] +git-tree-sha1 = "9ae9be75ad8ad9d26395bf625dea9beac6d519f1" +uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc" +version = "0.6.0" + +[[deps.Compat]] +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.16.0" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.1.1+0" + +[[deps.Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + +[[deps.DataAPI]] +git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.16.0" + +[[deps.DataFrames]] +deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] +git-tree-sha1 = "fb61b4812c49343d7ef0b533ba982c46021938a6" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "1.7.0" + +[[deps.DataFramesMeta]] +deps = ["Chain", "DataFrames", "MacroTools", "OrderedCollections", "Reexport", "TableMetadataTools"] +git-tree-sha1 = "21a4335f249f8b5f311d00d5e62938b50ccace4e" +uuid = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" +version = "0.15.4" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.20" + +[[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" +version = "1.11.0" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" +version = "1.11.0" + +[[deps.InlineStrings]] +git-tree-sha1 = "6a9fde685a7ac1eb3495f8e812c5a7c3711c2d5e" +uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" +version = "1.4.3" + + [deps.InlineStrings.extensions] + ArrowTypesExt = "ArrowTypes" + ParsersExt = "Parsers" + + [deps.InlineStrings.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" + Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +version = "1.11.0" + +[[deps.InvertedIndices]] +git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe" +uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +version = "1.3.1" + +[[deps.IteratorInterfaceExtensions]] +git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" +uuid = "82899510-4779-5014-852e-03e436cf321d" +version = "1.0.0" + +[[deps.LaTeXStrings]] +git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.4.0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +version = "1.11.0" + +[[deps.MacroTools]] +git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.15" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +version = "1.11.0" + +[[deps.Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.2.0" + +[[deps.OffsetArrays]] +git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.15.0" + + [deps.OffsetArrays.extensions] + OffsetArraysAdaptExt = "Adapt" + + [deps.OffsetArrays.weakdeps] + Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.27+1" + +[[deps.OrderedCollections]] +git-tree-sha1 = "cc4054e898b852042d7b503313f7ad03de99c3dd" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.0" + +[[deps.PooledArrays]] +deps = ["DataAPI", "Future"] +git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3" +uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" +version = "1.4.3" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.4.3" + +[[deps.PrettyTables]] +deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "1101cd475833706e4d0e7b122218257178f48f34" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.4.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SentinelArrays]] +deps = ["Dates", "Random"] +git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0" +uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" +version = "1.4.8" + +[[deps.SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.2.1" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] +git-tree-sha1 = "0feb6b9031bd5c51f9072393eb5ab3efd31bf9e4" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.9.13" + + [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" + StaticArraysStatisticsExt = "Statistics" + + [deps.StaticArrays.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.3" + +[[deps.Statistics]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0" +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.11.1" + + [deps.Statistics.extensions] + SparseArraysExt = ["SparseArrays"] + + [deps.Statistics.weakdeps] + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.StringManipulation]] +deps = ["PrecompileTools"] +git-tree-sha1 = "725421ae8e530ec29bcbdddbe91ff8053421d023" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.4.1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.TableMetadataTools]] +deps = ["DataAPI", "Dates", "TOML", "Tables", "Unitful"] +git-tree-sha1 = "c0405d3f8189bb9a9755e429c6ea2138fca7e31f" +uuid = "9ce81f87-eacc-4366-bf80-b621a3098ee2" +version = "0.1.0" + +[[deps.TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.1" + +[[deps.Tables]] +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] +git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "1.12.0" + +[[deps.Trapz]] +git-tree-sha1 = "79eb0ed763084a3e7de81fe1838379ac6a23b6a0" +uuid = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" +version = "2.0.3" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +version = "1.11.0" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" + +[[deps.Unitful]] +deps = ["Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "c0667a8e676c53d390a09dc6870b3d8d6650e2bf" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.22.0" + + [deps.Unitful.extensions] + ConstructionBaseUnitfulExt = "ConstructionBase" + InverseFunctionsUnitfulExt = "InverseFunctions" + + [deps.Unitful.weakdeps] + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.11.0+0" diff --git a/Project.toml b/Project.toml index b0873dd..4368344 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,25 @@ name = "CEF" uuid = "a0ebe571-4eea-4bd0-90ac-e9ca257b9c48" authors = ["J. A. Hernández"] -version = "1.0.0-DEV" +version = "0.1.0-DEV" + +[deps] +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" [compat] +DataFrames = "1.7.0" +DataFramesMeta = "0.15.4" +LinearAlgebra = "1.11.0" +OffsetArrays = "1.15.0" +StaticArrays = "1.9.13" +Statistics = "1.11.1" +Trapz = "2.0.3" julia = "1.7" [extras] diff --git a/README.md b/README.md index 07045a4..bdb2cf1 100644 --- a/README.md +++ b/README.md @@ -1,3 +1 @@ -# CEF - -[![Build Status](https://github.com/abehersan/CEF.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/abehersan/CEF.jl/actions/workflows/CI.yml?query=branch%3Amain) +# CEF.jl diff --git a/examples/hex_yb.jl b/examples/hex_yb.jl new file mode 100644 index 0000000..14408db --- /dev/null +++ b/examples/hex_yb.jl @@ -0,0 +1,62 @@ +""" +Reproduction of results from the paper by Rotter and Bauer +""" + + +using CEF +using DataFrames +using StatsPlots + + +function main() + + ion = single_ion("Yb3+") + bfactors = blm_dframe(Dict("B20"=>0.5622, "B40"=>1.6087e-5, + "B60"=>6.412e-7, "B66"=>-8.324e-6)) + + """ INS X-section """ + temps = [10.0, 50.0, 200.0] + ins_plot = plot(xlabel="Energy [meV]", ylabel="I(Q, E) [arb. units]") + for t in temps + calc_grid = DataFrame(T=t, Q=2.55, EN=range(-12,12,700)) + cef_neutronxsection_powder!(ion, bfactors, calc_grid) + @df calc_grid plot!(ins_plot, :EN, :I_CALC, label="T=$t K") + end + + """ Magnetic moment """ + mag_plot = plot(xlabel="Magnetic Field [Tesla]", ylabel="Magnetic Moment (muB)") + Bs = 0:0.5:12 + calc_grid = DataFrame(T=1.5, Bx=0.0, By=0.0, Bz=Bs) + cef_magneticmoment_crystal!(ion, bfactors, calc_grid, units=:ATOMIC) + @df calc_grid plot!(mag_plot, :Bz, :Mp_CALC, label="B parallel z") + calc_grid = DataFrame(T=1.5, Bx=Bs, By=0.0, Bz=0.0) + cef_magneticmoment_crystal!(ion, bfactors, calc_grid, units=:ATOMIC) + @df calc_grid plot!(mag_plot, :Bx, :Mp_CALC, label="B parallel x") + ylims!(mag_plot, 0, 3.5) + + """ Static susceptibility """ + chi_plot = plot(xlabel="Temperature [K]", ylabel="1/chi (emu/mol)") + calc_grid_z = DataFrame(T=1.5:0.5:300, Bx=0.0, By=0.0, Bz=0.01) + cef_susceptibility_crystal!(ion, bfactors, calc_grid_z, units=:CGS) + @df calc_grid_z plot!(chi_plot, :T, 1 ./ :CHIz_CALC, label="B parallel z") + calc_grid_x = DataFrame(T=1.5:0.5:300, Bx=0.01, By=0.0, Bz=0.0) + cef_susceptibility_crystal!(ion, bfactors, calc_grid_x, units=:CGS) + @df calc_grid_x plot!(chi_plot, :T, 1 ./ :CHIx_CALC, label="B parallel x") + calc_grid_powd = DataFrame(T=1.5:0.5:300) + cef_susceptibility_powder!(ion, bfactors, calc_grid_powd, units=:CGS) + @df calc_grid_powd plot!(chi_plot, :T, 1 ./ :CHI_CALC, label="Powder") + + """ Specific heat capacity and magnetic entropy """ + s_plot = plot(xlabel="Temperature [K]", ylabel="S, HC (J/mol/K)") + calc_grid = DataFrame(T=0.5:0.5:300, Bx=0.0, By=0.0, Bz=0.0) + cef_entropy!(ion, bfactors, calc_grid) + @df calc_grid plot!(s_plot, :T, :HC_CALC, label="HC") + @df calc_grid plot!(s_plot, :T, :SM_CALC, label="SM") + + plt = plot(ins_plot, mag_plot, chi_plot, s_plot, layout=(2, 2), size=(1080, 720)) + display(plt) + return nothing +end + + +main() \ No newline at end of file diff --git a/examples/re_neutronformfactor.jl b/examples/re_neutronformfactor.jl new file mode 100644 index 0000000..63967b6 --- /dev/null +++ b/examples/re_neutronformfactor.jl @@ -0,0 +1,39 @@ +""" +Plot the form factor for neutrons for all trivalent rare-earth ions +The x-axis is in Angstrom +The y-axis is the dipolar form factor squared calculated +""" + + +using CEF +using DataFrames +using StatsPlots + + +function main() + + mag_ions = [ + "Ce3+", "Pr3+", "Nd3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+" + ] + + plts = [] + c = :steelblue + QS = 0:0.01:20 + + for i in eachindex(mag_ions) + ion = single_ion(mag_ions[i]) + lbl = "Ion: $(ion.ion), J: $(ion.J)" + FSQUARED = abs2.(broadcast(q->CEF.dipolar_formfactor(ion, q), QS)) + plt = plot(xlims=(0, 15))#, ylims=(0, 1.5)) + plot!(QS, FSQUARED, label=lbl, c=c) + push!(plts, plt) + end + + plt = plot(plts..., layout=(5, 2), size=(1080, 720)) + display(plt) + + return nothing +end + + +main() \ No newline at end of file diff --git a/src/CEF.jl b/src/CEF.jl index fbb0c88..2d6cb2c 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -1,5 +1,68 @@ module CEF -# Write your package code here. +using DataFrames +using DataFramesMeta +using LinearAlgebra +using OffsetArrays +using StaticArrays +using Statistics +using Trapz -end +const PREC::Float64 = 1.0e-7 # for degeneracy calculations +const SDIG::Int64 = 7 # for numerical cutoffs + +# static vector/array definitions +const VEC3 = SVector{3, Float64} +const MAT3 = SMatrix{3, 3, Float64, 9} +const VEC{N} = SVector{N, Float64} +const MVEC{N} = MVector{N, Float64} +const CVEC{N} = SVector{N, ComplexF64} +const CMAT{N} = SMatrix{N, N, ComplexF64} +const HERMITIANC64 = Hermitian{ComplexF64, Matrix{ComplexF64}} + +# physical constants +const meV_per_K = 0.086173332621451774 # divide E in meV by meV_per_K to get E in K +const mu0 = 1.25663706212e-6 # [N/A] +const muB = 0.057883738013331 # [meV/T] +const kB = 0.08617333262 # [meV/K] +const NA = 6.02214076e23 # Avogadro constant, [1/mol] +const Rg = 8.314462618 # Ideal gas constant, [J/mol/K] + +include("./single_ion.jl") +export single_ion +export custom_ion +export mag_ion +export spin_operators +export re_hundsrules + +include("./cef_utils.jl") +export blm_dframe +export get_alm! + +include("./cef_system.jl") +export cef_diagonalization +export cef_hamiltonian +export stevens_EO, stevens_O + +include("./thermodynamical_quantities.jl") +export partition_function +export population_factor +export thermal_average + +include("./cef_entropy.jl") +export cef_entropy! +export cef_entropy_speclevels! + +include("./cef_magnetization.jl") +export cef_magneticmoment_crystal! +export cef_magneticmoment_powder! + +include("./cef_susceptibility.jl") +export cef_susceptibility_crystal! +export cef_susceptibility_powder! + +include("./cef_neutronxsection.jl") +export cef_neutronxsection_crystal! +export cef_neutronxsection_powder! + +end \ No newline at end of file diff --git a/src/cef_entropy.jl b/src/cef_entropy.jl new file mode 100644 index 0000000..a26e047 --- /dev/null +++ b/src/cef_entropy.jl @@ -0,0 +1,53 @@ +function hc_units(units::Symbol) + if isequal(units, :SI) + return Rg + else + @error "Units not supported in heat capacity calculations. Use SI, (J/K/mol)" + end +end + + +function mag_heatcap(Ep::Vector{Float64}, T::Real)::Float64 + np = population_factor(Ep, T) + heatcap = sum( ( Ep ./ (kB*T) ) .^2 .* np) - sum( Ep ./ (kB*T) .* np )^2 + return heatcap +end + + +function mag_entropy(HC::Vector{Float64}, T::Vector{Float64})::Vector{Float64} + S = similar(HC) + @views @inbounds for i in eachindex(T) + S[i] = trapz(T[1:i], HC[1:i] ./ T[1:i]) + end + return S +end + + +function cef_entropy!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; units::Symbol=:SI, method::Symbol=:EO)::Nothing + convfac = hc_units(units) + @eachrow! dfcalc begin + @newcol :HC_CALC::Vector{Float64} + B=[:Bx,:By,:Bz] + E=eigvals(cef_hamiltonian(ion,cefparams;B=B,method=method)) + E .-= minimum(E) + :HC_CALC=round(mag_heatcap(E,:T)*convfac,digits=SDIG) + end + dfcalc[:,:SM_CALC]=round.(mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]),digits=SDIG) + return nothing +end + + +function cef_entropy_speclevels!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; levels::UnitRange=1:4, units::Symbol=:SI, method::Symbol=:EO)::Nothing + # only levels specified contribute (2J+1 levels total) + convfac = hc_units(units) + @eachrow! dfcalc begin + @newcol :HC_CALC::Vector{Float64} + B=[:Bx,:By,:Bz] + E=eigvals(cef_hamiltonian(ion,cefparams;B=B,method=method)) + E .-= minimum(E) + E=E[levels] + :HC_CALC=round(mag_heatcap(E,:T)*convfac,digits=SDIG) + end + dfcalc[:,:SM_CALC]=round.(mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]),digits=SDIG) + return nothing +end \ No newline at end of file diff --git a/src/cef_magnetization.jl b/src/cef_magnetization.jl new file mode 100644 index 0000000..45558a5 --- /dev/null +++ b/src/cef_magnetization.jl @@ -0,0 +1,78 @@ +function mag_units(units::Symbol)::Float64 + if isequal(units, :SI) + return 5.5849397 # NA * muB ( J/T/mol ) + elseif isequal(units, :CGS) + return 5.5849397*1000.0 # NA * muB ( emu/mol ) + elseif isequal(units, :ATOMIC) + return 1.0 # units of Bohr magneton per mol + else + @error "Units $units not understood. Use one of either :SI, :CGS or :ATOMIC" + end +end + + +function calc_magmom(g::MVEC{3}, spinops::Vector{Matrix{ComplexF64}}, Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, T::Real, mode::Function)::Float64 + spin_expval = zeros(Float64, 3) + @inbounds for i in eachindex(spin_expval) + if iszero(spinops[i]) + continue + end + spin_expval[i] = thermal_average(Ep=Ep, Vp=Vp, op=spinops[i], T=T, mode=mode) + end + return g[1]*spin_expval[1]+g[2]*spin_expval[2]+g[3]*spin_expval[3] +end + + +function cef_magneticmoment_crystal!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) + unit_factor = mag_units(units) + spinops = [ion.Jx,ion.Jy,ion.Jz] + @eachrow! dfcalc begin + @newcol :Mp_CALC::Vector{Float64} + @newcol :Mx_CALC::Vector{Float64} + @newcol :My_CALC::Vector{Float64} + @newcol :Mz_CALC::Vector{Float64} + extfield = [:Bx, :By, :Bz] + E, V = eigen(cef_hamiltonian(ion,cefparams;B=extfield,method=method)) + E .-= minimum(E) + if iszero(extfield) + spin_proj = spinops + else + spin_proj = spinops .* normalize(extfield) + end + spin_expval = zeros(Float64, 3) + @inbounds for i in eachindex(spin_expval) + if iszero(spinops[i]) + continue + end + spin_expval[i] = thermal_average(Ep=E,Vp=V,op=spinops[i],T=:T,mode=mode) + end + :Mp_CALC=round(norm(ion.g .* spin_expval)*unit_factor,digits=SDIG) + :Mx_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[1],T=:T,mode=mode)*ion.g[1]*unit_factor,digits=SDIG) + :My_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[2],T=:T,mode=mode)*ion.g[2]*unit_factor,digits=SDIG) + :Mz_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[3],T=:T,mode=mode)*ion.g[3]*unit_factor,digits=SDIG) + end + return nothing +end + + +function cef_magneticmoment_powder!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) + unit_factor = mag_units(units) + spinops = [ion.Jx,ion.Jy,ion.Jz] + @eachrow! dfcalc begin + @newcol :M_CALC::Vector{Float64} + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[:B,0.0,0.0],method=method)) + E .-= minimum(E) + MX = calc_magmom(ion.g,spinops .* [1.0, 0.0, 0.0],E,V,:T,mode) + + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[0.0,:B,0.0],method=method)) + E .-= minimum(E) + MY = calc_magmom(ion.g,spinops .* [0.0, 1.0, 0.0],E,V,:T,mode) + + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[0.0,0.0,:B],method=method)) + E .-= minimum(E) + MZ = calc_magmom(ion.g,spinops .* [0.0, 0.0, 1.0],E,V,:T,mode) + + :M_CALC=round(((MX + MY + MZ) / 3.0) * unit_factor,digits=SDIG) + end + return nothing +end \ No newline at end of file diff --git a/src/cef_neutronxsection.jl b/src/cef_neutronxsection.jl new file mode 100644 index 0000000..7214412 --- /dev/null +++ b/src/cef_neutronxsection.jl @@ -0,0 +1,144 @@ +function dipolar_formfactor(ion::mag_ion, Q::Real)::Float64 + A_j0, a_j0, B_j0, b_j0, C_j0, c_j0, D_j0 = ion.ff_coeff_j0 + A_j2, a_j2, B_j2, b_j2, C_j2, c_j2, D_j2 = ion.ff_coeff_j2 + s = Q / 4pi + ff_j0 = A_j0 * exp(-a_j0*s^2) + B_j0 * exp(-b_j0*s^2) + C_j0 * exp(-c_j0*s^2) + D_j0 + ff_j2 = A_j2*s^2 * exp(-a_j2*s^2) + B_j2*s^2 * exp(-b_j2*s^2) + C_j2*s^2 * exp(-c_j2*s^2) + D_j2*s^2 + return ff_j0 + ion.C2 * ff_j2 +end + + +function calc_polmatrix(Qcart::Vector{Float64})::Matrix{Float64} + polmat = Matrix{Float64}(undef, (3, 3)) + Q = normalize(Qcart) + for I in CartesianIndices(polmat) + a, b = Tuple(I) + polmat[I] = (isequal(a, b)*1.0 - Q[a]*Q[b]) + end + return polmat +end + + +function calc_transitions(ion::mag_ion, i::Int64, Vp::Matrix{ComplexF64})::Vector{VEC{9}} + SIGMAS = VEC{9}[] + @views @inbounds for j in 1:size(Vp, 1) + sxx = real( dot( Vp[:, i], ion.Jx, Vp[:, j] ) * dot( Vp[:, j], ion.Jx, Vp[:, i] ) ) * (ion.g[1]*ion.g[1]) + sxy = real( dot( Vp[:, i], ion.Jy, Vp[:, j] ) * dot( Vp[:, j], ion.Jx, Vp[:, i] ) ) * (ion.g[1]*ion.g[2]) + sxz = real( dot( Vp[:, i], ion.Jz, Vp[:, j] ) * dot( Vp[:, j], ion.Jx, Vp[:, i] ) ) * (ion.g[1]*ion.g[3]) + syx = real( dot( Vp[:, i], ion.Jx, Vp[:, j] ) * dot( Vp[:, j], ion.Jy, Vp[:, i] ) ) * (ion.g[2]*ion.g[1]) + syy = real( dot( Vp[:, i], ion.Jy, Vp[:, j] ) * dot( Vp[:, j], ion.Jy, Vp[:, i] ) ) * (ion.g[2]*ion.g[2]) + syz = real( dot( Vp[:, i], ion.Jz, Vp[:, j] ) * dot( Vp[:, j], ion.Jy, Vp[:, i] ) ) * (ion.g[2]*ion.g[3]) + szx = real( dot( Vp[:, i], ion.Jx, Vp[:, j] ) * dot( Vp[:, j], ion.Jz, Vp[:, i] ) ) * (ion.g[3]*ion.g[1]) + szy = real( dot( Vp[:, i], ion.Jy, Vp[:, j] ) * dot( Vp[:, j], ion.Jz, Vp[:, i] ) ) * (ion.g[3]*ion.g[2]) + szz = real( dot( Vp[:, i], ion.Jz, Vp[:, j] ) * dot( Vp[:, j], ion.Jz, Vp[:, i] ) ) * (ion.g[3]*ion.g[3]) + push!(SIGMAS, [sxx, sxy, sxz, syx, syy, syz, szx, szy, szz]) + end + return SIGMAS +end + + +function calc_neutronspectrum_xtal(ion::mag_ion, Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, Qcart::Vector{<:Real}, T::Real)::Vector{VEC{2}} + np = population_factor(Ep, T) + ffactor = dipolar_formfactor(ion, norm(Qcart)) + polfactors = reshape(calc_polmatrix(Qcart), 9) + NXS = VEC{2}[] + @inbounds for i in eachindex(Ep) + if isapprox(np[i], 0.0, atol=PREC) + continue + end + sigmas = calc_transitions(ion, i, Vp) + @inbounds for j in eachindex(Ep) + dE = Ep[j] - Ep[i] + NINT = dot(polfactors, sigmas[j])*np[i]*abs2(ffactor) + if dE < 0.0 # detailed balance + NINT *= exp( -abs(dE)/(kB*T) ) + end + push!(NXS, [dE, NINT]) + end + end + return NXS +end + + +function calc_neutronspectrum_powd(ion::mag_ion, Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, Qlen::Real, T::Real)::Vector{VEC{2}} + np = population_factor(Ep, T) + ffactor = dipolar_formfactor(ion, Qlen) + NXS = VEC{2}[] + @inbounds for i in eachindex(Ep) + if isapprox(np[i], 0.0, atol=PREC) + continue + end + sigmas = calc_transitions(ion, i, Vp) + @inbounds for j in eachindex(Ep) + dE = Ep[j] - Ep[i] + NINT = sum(sigmas[j])*np[i]*abs2(ffactor)*2.0/3.0 + if dE < 0.0 # detailed balance + NINT *= exp( -abs(dE)/(kB*T) ) + end + push!(NXS, [dE, NINT]) + end + end + return NXS +end + + +function simulate_Escan(NXS::Vector{VEC{2}}, Es::AbstractVector, R::Function=TAS_resfunc)::Vector{Float64} + Is = zeros(length(Es)) + @inbounds for i in eachindex(Es) + @inbounds for j in eachindex(NXS) + E, I = NXS[j] + Is[i] += I * R(Es[i], E) + end + end + return Is +end + + +function cef_neutronxsection_crystal!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; resfunc::Function=TAS_resfunc, method::Symbol=:EO)::Nothing + extfield = [mean(dfcalc.Bx), mean(dfcalc.By), mean(dfcalc.Bz)] + E, V = eigen(cef_hamiltonian(ion,cefparams,B=extfield,method=method)) + E .-= minimum(E) + T = mean(dfcalc.T) + Q = [mean(dfcalc.Qx), mean(dfcalc.Qy), mean(dfcalc.Qz)] + NINT = calc_neutronspectrum_xtal(ion,E,V,Q,T) + EN = dfcalc.EN + II = simulate_Escan(NINT, EN, resfunc) + dfcalc[!, :I_CALC] = II + return nothing +end + + +function cef_neutronxsection_powder!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; resfunc::Function=TAS_resfunc, method::Symbol=:EO)::Nothing + E, V = eigen(cef_hamiltonian(ion,cefparams;method=method)) + E .-= minimum(E) + T = mean(dfcalc.T) + Q = mean(dfcalc.Q) + NINT = calc_neutronspectrum_powd(ion,E,V,Q,T) + EN = dfcalc.EN + II = simulate_Escan(NINT, EN, resfunc) + dfcalc[!, :I_CALC] = II + return nothing +end + + +@doc raw""" + TAS_resfunc(E::Float64, Epeak::Float64, width::Function=x->0.2/(2*sqrt(2*log(2))))::Float64 + +Default resolution function. +This function evaluates a normalized Gaussian profile centered at `Epeak`. +The width of the peak is a function of `E`. +In this default, it is assumed constant and given by `width`. +""" +function TAS_resfunc(E::Float64, Epeak::Float64, width::Function=x->0.2/(2*sqrt(2*log(2))))::Float64 + return gaussian(x=E, mu=Epeak, A=1.0, sigma=width(E)) +end + + +function gaussian(; x::Real, A::Real, mu::Real, sigma::Real)::Float64 + return (A/(sqrt(2pi)*sigma))*exp(-(x-mu)^2/(2*sigma^2)) +end + + +function lorentz(; x::Real, A::Real, mu::Real, gamma::Real)::Float64 + return (A/pi)*(gamma/((x-mu)^2+gamma^2)) +end \ No newline at end of file diff --git a/src/cef_susceptibility.jl b/src/cef_susceptibility.jl new file mode 100644 index 0000000..e4f186c --- /dev/null +++ b/src/cef_susceptibility.jl @@ -0,0 +1,91 @@ +function chi_units(units::Symbol)::Float64 + if isequal(units, :SI) + return 4.062426*1e-7 # N_A * muB(erg/G) * muB(meV/G) + elseif isequal(units, :CGS) + return 0.03232776 # N_A * muB(J/T) * muB(meV/T) * mu0 + elseif isequal(units, :ATOMIC) + return 1.0 # muB / magnetic ion + else + @error "Units $units not understood. Use one of either :SI, :CGS or :ATOMIC" + end +end + + +function calc_chialphaalpha(; op_alpha::Matrix{ComplexF64}, Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, T::Real, mode::Function=real)::Float64 + chi_alphaalpha::Float64 = 0.0 + np = population_factor(Ep, T) + @views @inbounds for p in eachindex(Ep) + ep = Ep[p] + @views @inbounds for pp in eachindex(Ep) + epp = Ep[pp] + if isapprox(ep, epp; atol=PREC) + if isapprox(np[p], 0.0, atol=PREC) + continue + else + m_element_alpha = dot(Vp[:,p], op_alpha, Vp[:,pp]) + m_element = m_element_alpha*conj(m_element_alpha) + chi_alphaalpha += m_element*np[p]/(kB*T) + end + else + if isapprox(np[p], 0.0, atol=PREC) && isapprox(np[pp], 0.0, atol=PREC) + continue + else + m_element_alpha = dot(Vp[:,p], op_alpha, Vp[:,pp]) + m_element = m_element_alpha*conj(m_element_alpha) + pop_diff = np[p] - np[pp] + chi_alphaalpha += (m_element*pop_diff)/(epp-ep) + end + end + end + end + t_avg_alpha = thermal_average(Ep=Ep,Vp=Vp,op=op_alpha,T=T,mode=mode) + chi_alphaalpha -= (t_avg_alpha^2)/(kB*T) + return chi_alphaalpha +end + + +function calc_chi(g::MVEC{3}; spin_proj::Vector{Matrix{ComplexF64}}, Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, T::Real, mode::Function)::Float64 + chi_vec = zeros(Float64, 3) + @inbounds for i in eachindex(spin_proj) + if iszero(spin_proj[i]) + continue + else + chi_vec[i]=calc_chialphaalpha(op_alpha=spin_proj[i],Ep=Ep,Vp=Vp,T=T,mode=mode) + end + end + return norm((g .^2) .* chi_vec) +end + + +function cef_susceptibility_crystal!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; units::Symbol=:CGS, method::Symbol=:EO, mode::Function=real)::Nothing + unit_factor = chi_units(units) + spin_ops = [ion.Jx,ion.Jy,ion.Jz] + extfield = [mean(dfcalc.Bx), mean(dfcalc.By), mean(dfcalc.Bz)] + E, V = eigen(cef_hamiltonian(ion,cefparams;B=extfield,method=method)) + E .-= minimum(E) + @eachrow! dfcalc begin + @newcol :CHIx_CALC::Vector{Float64} + @newcol :CHIy_CALC::Vector{Float64} + @newcol :CHIz_CALC::Vector{Float64} + :CHIx_CALC=round(calc_chi(ion.g,spin_proj=spin_ops.*[1.0,0.0,0.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor,digits=SDIG) + :CHIy_CALC=round(calc_chi(ion.g,spin_proj=spin_ops.*[0.0,1.0,0.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor,digits=SDIG) + :CHIz_CALC=round(calc_chi(ion.g,spin_proj=spin_ops.*[0.0,0.0,1.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor,digits=SDIG) + end + return nothing +end + + +function cef_susceptibility_powder!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; units::Symbol=:CGS, method::Symbol=:EO, mode::Function=real)::Nothing + unit_factor = chi_units(units) + spin_ops = [ion.Jx,ion.Jy,ion.Jz] + E, V = eigen(cef_hamiltonian(ion,cefparams;method=method)) + E .-= minimum(E) + @eachrow! dfcalc begin + @newcol :CHI_CALC::Vector{Float64} + chix = calc_chi(ion.g,spin_proj=spin_ops.*[1.0,0.0,0.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor + chiy = calc_chi(ion.g,spin_proj=spin_ops.*[0.0,1.0,0.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor + chiz = calc_chi(ion.g,spin_proj=spin_ops.*[0.0,0.0,1.0],Ep=E,Vp=V,T=:T,mode=mode)*unit_factor + :CHI_CALC=round(((chix+chiy+chiz)/3),digits=SDIG) + end + return nothing +end \ No newline at end of file diff --git a/src/cef_system.jl b/src/cef_system.jl new file mode 100644 index 0000000..3a73ac8 --- /dev/null +++ b/src/cef_system.jl @@ -0,0 +1,259 @@ +function cef_diagonalization(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::Nothing + cef_matrix = cef_hamiltonian(ion,cefparams;B=B,method=method) + @assert ishermitian(cef_matrix) + E = eigvals(cef_matrix) + E .-= minimum(E) + printstyled("CEF matrix diagonalization results.\n\n", color=:underline, bold=true) + display(ion) + println("Diagonal g-tensor [gxx, gyy, gzz]: $(ion.g).\n") + println("External magnetic field in (Tesla) [Bx, By, Bz]: $B\n") + println("CEF energy levels in (meV) and in (K):") + for i in eachindex(E) + Emev = round(E[i], digits=SDIG) + EK = round(E[i]/meV_per_K, digits=SDIG) + Es = Emev, EK + println(join(Es, ",\t")) + end + return nothing +end + + +function cef_hamiltonian(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::HERMITIANC64 + if iszero(B) + return H_cef(ion, cefparams, method) + else + return H_cef(ion, cefparams, method) + H_zeeman(ion, B) + end +end + + +function cef_hamiltonian(ion::mag_ion, D::Real, E::Real; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::HERMITIANC64 + m_dim = Int(2*ion.J+1) + h_cef = zeros(ComplexF64, (m_dim, m_dim)) + h_cef += D * ion.Jz^2 + h_cef += E * (ion.Jp^2 + ion.Jm^2)/2 + if iszero(B) + return Hermitian(h_cef) + else + return Hermitian(h_cef) + H_zeeman(ion, B) + end +end + + +function H_zeeman(ion::mag_ion, extfield::Vector{<:Real})::HERMITIANC64 + Bx, By, Bz = extfield + gxx, gyy, gzz = ion.g + gxxBxJx = (gxx * Bx) * ion.Jx + gyyByJy = (gyy * By) * ion.Jy + gzzBzJz = (gzz * Bz) * ion.Jz + zeeman_matrix = (-1.0 * muB) * (gxxBxJx .+ gyyByJy .+ gzzBzJz) + return Hermitian(zeeman_matrix) +end + + +function H_cef(ion::mag_ion, cefparams::DataFrame, method::Symbol)::HERMITIANC64 + m_dim = Int(2*ion.J+1) + cef_matrix = zeros(ComplexF64, (m_dim, m_dim)) + if isequal(method,:EO) + mode=stevens_EO + elseif isequal(method,:O) + mode=stevens_O + else + @error "Method: $(method) not one of `:EO` or `:O`." + end + @eachrow! cefparams begin + cef_matrix += :B * mode(ion, :l, :m) + end + return Hermitian(cef_matrix) +end + + +function ryabov_clm(l::Int, m::Int)::Float64 + lmax = 13 + if l > lmax + @error "Invalid l, l= 0 + Op = clm/2 * (T + adjoint(T)) + else + Op = clm/2im * (T - adjoint(T)) + end + return Hermitian(Op) +end + + +function stevens_O(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 + J=ion.J + m_dim = Int(2*J+1) + Jz = ion.Jz + Jp = ion.Jp + Jm = ion.Jm + X = Diagonal(fill(J*(J+1), m_dim)) + if l == 1 + if m == 1 + O11 = 1.0/2.0 * (Jp + Jm) + return Hermitian(O11) + elseif m == -1 + O1m1 = -1.0im/2.0 * (Jp - Jm) + return Hermitian(O1m1) + elseif m == 0 + O10 = Jz + return Hermitian(O10) + else + err_message = + "Given values of l=$l and m=$m not supported.\n"* + "Only l=[2, 4, 6] and m = -l:1:l are currently supported." + @error err_message + end + elseif l == 2 + if m == -2 + O2m2 = -1.0im/2.0 * (Jp^2 - Jm^2) + return Hermitian(O2m2) + elseif m == -1 + O2m1 = -1.0im/4.0 * (Jz*(Jp - Jm) + (Jp - Jm)*Jz) + return Hermitian(O2m1) + elseif m == 0 + O20 = 3.0*Jz^2 - X + return Hermitian(O20) + elseif m == 1 + O21 = 1.0/4.0 * (Jz*(Jp + Jm) + (Jp + Jm)*Jz) + return Hermitian(O21) + elseif m == 2 + O22 = 1.0/2.0 * (Jp^2 + Jm^2) + return Hermitian(O22) + else + err_message = + "Given values of l=$l and m=$m not supported.\n"* + "Only l=[2, 4, 6] and m = -l:1:l are currently supported." + @error err_message + end + elseif l == 4 + if m == -4 + O4m4 = -1.0im/2.0 * (Jp^4 - Jm^4) + return Hermitian(O4m4) + elseif m == -3 + O4m3 = -1.0im/4.0 * ((Jp^3 - Jm^3)*Jz + Jz*(Jp^3 - Jm^3)) + return Hermitian(O4m3) + elseif m == -2 + O4m2 = -1.0im/4.0 * ((Jp^2 - Jm^2)*(7.0*Jz^2 - X - 5.0I) + (7.0*Jz^2 - X - 5.0I)*(Jp^2 - Jm^2)) + return Hermitian(O4m2) + elseif m == -1 + O4m1 = -1.0im/4.0 * ((Jp - Jm)*(7.0*Jz^3 - (3.0*X + 1.0I)*Jz) + (7.0*Jz^3 - (3.0*X + 1.0I)*Jz)*(Jp - Jm)) + return Hermitian(O4m1) + elseif m == 0 + O40 = 35.0*Jz^4 - (30.0*X - 25.0I)*Jz^2 + 3.0*X^2 - 6.0*X + return Hermitian(O40) + elseif m == 1 + O41 = 1.0/4.0 * ((Jp + Jm)*(7.0*Jz^3 - (3.0*X + 1.0I)*Jz) + (7.0*Jz^3 - (3.0*X + 1.0I)*Jz)*(Jp + Jm)) + return Hermitian(O41) + elseif m == 2 + O42 = 1.0/4.0 * ((Jp^2 + Jm^2)*(7.0*Jz^2 - X - 5.0I) + (7.0*Jz^2 - X - 5.0I)*(Jp^2 + Jm^2)) + return Hermitian(O42) + elseif m == 3 + O43 = 1.0/4.0 * (Jz*(Jp^3 + Jm^3) + (Jp^3 + Jm^3)*Jz) + return Hermitian(O43) + elseif m == 4 + O44 = 1.0/2.0 * (Jp^4 + Jm^4) + return Hermitian(O44) + else + err_message = + "Given values of l=$l and m=$m not supported.\n"* + "Only l=[2, 4, 6] and m = -l:1:l are currently supported." + @error err_message + end + elseif l == 6 + if m == -6 + O6m6 = -1.0im/2.0 * (Jp^6 - Jm^6) + return Hermitian(O6m6) + elseif m == -5 + O6m5 = -1.0im/4.0 * ((Jp^5 - Jm^5)*Jz + Jz*(Jp^5 - Jm^5)) + return Hermitian(O6m5) + elseif m == -4 + O6m4 = -1.0im/4.0 * ((Jp^4 - Jm^4)*(11.0*Jz^2 - X - 38.0I) + (11.0*Jz^2 - X - 38.0I)*(Jp^4 - Jm^4)) + return Hermitian(O6m4) + elseif m == -3 + O6m3 = -1.0im/4.0 * ((Jp^3 - Jm^3)*(11.0*Jz^3 - (3.0*X + 59.0I)*Jz) + (11.0*Jz^3 - (3.0*X + 59.0I)*Jz)*(Jp^3 - Jm^3)) + return Hermitian(O6m3) + elseif m == -2 + O6m2 = -1.0im/4.0 * ((Jp^2 - Jm^2)*(33.0*Jz^4 - (18.0*X + 123.0I)*Jz^2 + X^2 + 10.0*X + 102I) + (33.0*Jz^4 - (18.0*X + 123.0I)*Jz^2 + X^2 + 10.0*X + 102.0I)*(Jp^2 - Jm^2)) + return Hermitian(O6m2) + elseif m == -1 + O6m1 = -1.0im/4.0 * ((Jp - Jm)*(33.0*Jz^5 - (30.0*X - 15.0I)*Jz^3 + (5.0*X^2 - 10.0*X + 12.0I)*Jz) + (33.0*Jz^5 - (30.0*X - 15.0I)*Jz^3 + (5.0*X^2 - 10.0*X + 12.0I)*Jz)*(Jp - Jm)) + return Hermitian(O6m1) + elseif m == 0 + O60 = 231.0*Jz^6 - (315.0*X - 735.0I)*Jz^4 + (105.0*X^2 - 525.0*X + 294.0I)*Jz^2 - 5.0*X^3 + 40.0*X^2 - 60.0*X + return Hermitian(O60) + elseif m == 1 + O61 = 1.0/4.0 * ((Jp + Jm)*(33.0*Jz^5 - (30.0X - 15.0I)*Jz^3 + (5.0*X^2 - 10.0*X + 12.0I)*Jz) + (33.0*Jz^5 - (30.0*X - 15.0I)*Jz^3 + (5.0*X^2 - 10.0*X + 12.0I)*Jz)*(Jp + Jm)) + return Hermitian(O61) + elseif m == 2 + O62 = 1.0/4.0 *((Jp^2 + Jm^2)*(33.0*Jz^4 - (18*X + 123.0I)*Jz^2 + X^2 + 10.0*X + 102.0I) + (33.0*Jz^4 - (18.0*X + 123.0I)*Jz^2 + X^2 + 10.0*X + 102.0I)*(Jp^2 + Jm^2)) + return Hermitian(O62) + elseif m == 3 + O63 = 1.0/4.0 * ((Jp^3 + Jm^3)*(11.0*Jz^3 - (3.0*X + 59.0I)*Jz) + (11.0*Jz^3 - (3.0*X + 59.0I)*Jz)*(Jp^3 + Jm^3)) + return Hermitian(O63) + elseif m == 4 + O64 = 1.0/4.0 * ((Jp^4 + Jm^4)*(11.0*Jz^2 - X - 38.0I) + (11.0*Jz^2 - X - 38.0I)*(Jp^4 + Jm^4)) + return Hermitian(O64) + elseif m == 5 + O65 = 1.0/4.0 * ((Jp^5 + Jm^5)*Jz + Jz*(Jp^5 + Jm^5)) + return Hermitian(O65) + elseif m == 6 + O66 = 1.0/2.0 * (Jp^6 + Jm^6) + return Hermitian(O66) + else + err_message = + "Given values of l=$l and m=$m not supported.\n"* + "Only l=[2, 4, 6] and m = -l:1:l are currently supported." + @error err_message + end + else + err_message = + "Given values of l=$l and m=$m not supported.\n"* + "Only l=[2, 4, 6] and m = -l:1:l are currently supported." + @error err_message + end +end \ No newline at end of file diff --git a/src/cef_utils.jl b/src/cef_utils.jl new file mode 100644 index 0000000..3970fc3 --- /dev/null +++ b/src/cef_utils.jl @@ -0,0 +1,58 @@ +function blm_dframe(blm_dict::Dict{String, <:Real})::DataFrame + l = Int[] + m = Int[] + B = Float64[] + for k in keys(blm_dict) + ll, mm = parse_blm(k) + BB = blm_dict[k] + push!(l, ll) + push!(m, mm) + push!(B, BB) + end + return DataFrame("B"=>B, "l"=>l, "m"=>m) +end + + +function parse_blm(b::String)::Tuple{Int, Int} + if b[3] == 'm' + l = parse(Int, b[2]) + m = -1*parse(Int, b[end]) + else + l = parse(Int, b[2]) + m = parse(Int, b[end]) + end + return (l, m) +end + + +@doc raw""" + get_alm!(bfactors::DataFrame, single_ion::mag_ion)::DataFrame + +Factorization of the Stevens B_lm parameters defined as +B_lm = A_lm * * theta_l, +where is the expectation value of the radial wavefunction of the +4f electron density (tabulated) and theta_l are the Stevens geometrical factors. +""" +function get_alm!(bfactors::DataFrame, single_ion::mag_ion) + alpha, beta, gamma = single_ion.stevens_factors + r2, r4, r6 = single_ion.rad_wavefunction + alm = zeros(Float64, nrow(bfactors)) + for (i, r) in enumerate(eachrow(bfactors)) + if r.l == 2 + alm[i] = r.B / (alpha * r2) + elseif r.l == 4 + alm[i] = r.B / (beta * r4) + elseif r.l == 6 + alm[i] = r.B / (gamma * r6) + else + alm[i] = r.B + err_message = + "Given Blm parameter has invalid l and/or m.\n"* + "l=$l and m=$m were parsed and are not supported.\n"* + "Writing unscaled Stevens parameter." + @error err_message + end + end + bfactors[!, :A] .= alm + return nothing +end \ No newline at end of file diff --git a/src/single_ion.jl b/src/single_ion.jl new file mode 100644 index 0000000..f5c2600 --- /dev/null +++ b/src/single_ion.jl @@ -0,0 +1,190 @@ +const re_hundsrules = Dict( + "Ce3+"=>[ + 5/2, 6/7, # J, gJ + -5.7143/1e2, 63.4921/1e4, 0, # , , + 0.3666, 0.3108, 0.5119, # Stevens geometrical factors + 0.2291, 18.18, 0.7897, 5.807, -0.0191, 0.0, 0.0, # for j0 + 2.1284, 8.9174, 1.1229, 2.8371, 0.01108, 0.0, 0.0 # for j2 + ], + "Pr3+"=>[ + 4/1, 4/5, + -2.101/1e2, -7.3462/1e4, 60.994/1e6, + 0.3350, 0.2614, 0.4030, + 0.0504, 24.9989, 0.2572, 12.0377, 0.7142, 5.0039, -0.0219, + 0.8734, 18.9876, 1.5594, 6.0872, 0.8142, 2.4150, 0.0111 + ], + "Nd3+"=>[ + 9/2, 8/11, + -0.6428/1e2, -2.9111/1e4, -37.9980/1e6, + 0.3120, 0.2282, 0.3300, + 0.0540, 25.0293, 0.3101, 12.1020, 0.6575, 4.7223, -0.0216, + 0.6751, 18.3421, 1.6272, 7.2600, 0.9644, 2.6016, 0.0150 + ], + "Pm3+"=>[ + 4/1, 3/5, + 0.7713/1e2, 4.0755/1e4, 60.7807/1e6, + 0.2899, 0.1991, 0.2755, + 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0 + ], + "Sm3+"=>[ + 5/2, 2/7, + 4.127/1e2, 25.0120/1e4, 0, + 0.2728, 0.1772, 0.2317, + 0.0288, 25.2068, 0.2973, 11.8311, 0.6954, 4.2117, -0.0213, + 0.4707, 18.4301, 1.4261, 7.0336, 0.9574, 2.4387, 0.0182 + ], + "Tb3+"=>[ + 6/1, 3/2, + -1.0101/1e2, 1.2244/1e4, -1.1212/1e6, + 0.2302, 0.1295, 0.1505, + 0.0177, 25.5095, 0.2921, 10.5769, 0.7133, 3.5122, -0.0231, + 0.2892, 18.4973, 1.1678, 6.7972, 0.9437, 2.2573, 0.0232 + ], + "Dy3+"=>[ + 15/2, 4/3, + -0.6349/1e2, -0.592/1e4, 1.035/1e6, + 0.2188, 0.1180, 0.1328, + 0.1157, 15.0732, 0.3270, 6.7991, 0.5821, 3.0202, -0.0249, + 0.2523, 18.5172, 1.0914, 6.7362, 0.9345, 2.2082, 0.0250 + ], + "Ho3+"=>[ + 8/1, 5/4, + -0.2222/1e2, -0.333/1e4, -1.2937/1e6, + 0.2085, 0.1081, 0.1181, + 0.0566, 18.3176, 0.3365, 7.6880, 0.6317, 2.9427, -0.0248, + 0.2188, 18.5157, 1.0240, 6.7070, 0.9251, 2.1614, 0.0268 + ], + "Er3+"=>[ + 15/2, 6/5, + 0.254/1e2, 0.444/1e4, 2.0699/1e6, + 0.1991, 0.0996, 0.1058, + 0.0586, 17.9802, 0.3540, 7.0964, 0.6126, 2.7482, -0.0251, + 0.1710, 18.5337, 0.9879, 6.6246, 0.9044, 2.1004, 0.0278 + ], + "Tm3+"=>[ + 6/1, 7/6, + 1.0101/1e2, 1.632/1e4, -5.606/1e6, + 0.1905, 0.0921, 0.0953, + 0.0581, 15.0922, 0.2787, 7.8015, 0.6854, 2.7931, -0.0224, + 0.1760, 18.5417, 0.9105, 6.5787, 0.8970, 2.0622, 0.0294 + ], + "Yb3+"=>[ + 7/2, 8/7, + 3.175/1e2, -17.3160/1e4, 148.0/1e6, + 0.1826, 0.0854, 0.0863, + 0.0416, 16.0949, 0.2849, 7.8341, 0.6961, 2.6725, -0.0229, + 0.1570, 18.5553, 0.8484, 6.5403, 0.8880, 2.0367, 0.0318 + ], +) + + +function single_ion(ion::String) + try + J, g, + alpha, beta, gamma, + r2, r4, r6, + A_j0, a_j0, B_j0, b_j0, C_j0, c_j0, D_j0, + A_j2, a_j2, B_j2, b_j2, C_j2, c_j2, D_j2 = re_hundsrules[ion] + C2=(2-g)/g + return mag_ion( + ion, J, + spin_operators(J, "x"), + spin_operators(J, "y"), + spin_operators(J, "+"), + spin_operators(J, "-"), + spin_operators(J, "z"), + [g, g, g], + [alpha, beta, gamma], + [r2, r4, r6], + C2, + [A_j0, a_j0, B_j0, b_j0, C_j0, c_j0, D_j0], + [A_j2, a_j2, B_j2, b_j2, C_j2, c_j2, D_j2], + ) + catch y + mag_ions = collect(keys(re_hundsrules)) + err_message = + "$y\n"* + "Given ion $ion not supported.\n"* + "Available magnetic ions: $mag_ions" + @error err_message + end +end + + +function custom_ion(S::Real, g::Real, C2::Real)::mag_ion + return mag_ion( + ion="S", J=S, + Jx=spin_operators(S,"x"), + Jy=spin_operators(S,"y"), + Jz=spin_operators(S,"z"), + Jp=spin_operators(S,"+"), + Jm=spin_operators(S,"-"), + g=[g,g,g], + stevens_factors=zeros(3), + rad_wavefunction=zeros(3), + C2=C2, + ff_coeff_j0=zeros(7), + ff_coeff_j2=zeros(7), + ) +end + + +Base.@kwdef mutable struct mag_ion + ion::String + J::Float64 + Jx::Matrix{ComplexF64} + Jy::Matrix{ComplexF64} + Jp::Matrix{ComplexF64} + Jm::Matrix{ComplexF64} + Jz::Matrix{ComplexF64} + g::MVEC{3} + stevens_factors::VEC{3} + rad_wavefunction::VEC{3} + C2::Float64 + ff_coeff_j0::VEC{7} + ff_coeff_j2::VEC{7} +end + + +function Base.show(io::IO, ::MIME"text/plain", ion::mag_ion) + printstyled(io, "Magnetic ion: $(ion.ion)\n") + if isequal(ion.ion, "S") + @warn "Custom ion! Geometric factors, radial wavefunction coefficients and form-factor coefficients are zero. Manually initialize them." + end + print(io, "Quantum number J: $(ion.J).\nHilbert space dimension: $(Int(2*ion.J+1)).") + return nothing +end + + +function spin_operators(J::Float64, a::String)::Matrix{ComplexF64} + mJ = -J:1:J + if isequal(a, "x") + jp_eigval = @. sqrt(J*(J+1)-mJ*(mJ+1)) + jm_eigval = @. sqrt(J*(J+1)-mJ*(mJ-1)) + Jp = diagm(1=>jp_eigval[1:end-1]) + Jm = diagm(-1=>jm_eigval[2:end]) + Jx = (Jp + Jm)/2.0 + return Jx + elseif isequal(a, "y") + jp_eigval = @. sqrt(J*(J+1)-mJ*(mJ+1)) + jm_eigval = @. sqrt(J*(J+1)-mJ*(mJ-1)) + Jp = diagm(1=>jp_eigval[1:end-1]) + Jm = diagm(-1=>jm_eigval[2:end]) + Jy = (Jp - Jm)/2.0im + return Jy + elseif isequal(a, "z") + Jz = diagm(mJ) + return Jz + elseif isequal(a, "+") + jp_eigval = @. sqrt(J*(J+1)-mJ*(mJ+1)) + Jp = diagm(1=>jp_eigval[1:end-1]) + return Jp + elseif isequal(a, "-") + jm_eigval = @. sqrt(J*(J+1)-mJ*(mJ-1)) + Jm = diagm(-1=>jm_eigval[2:end]) + return Jm + else + @error "String $(a) not understood. Choose one of either {x, y, z, +, -}" + end +end \ No newline at end of file diff --git a/src/thermodynamical_quantities.jl b/src/thermodynamical_quantities.jl new file mode 100644 index 0000000..f8ca7c7 --- /dev/null +++ b/src/thermodynamical_quantities.jl @@ -0,0 +1,42 @@ +function partition_function(Ep::Vector{Float64}, T::Real)::Float64 + if isapprox(T,0.0,atol=PREC) + return 1.0 + else + Z = 0.0 + @inbounds for i in eachindex(Ep) + Z += exp( -Ep[i]/(kB*T) ) + end + return Z + end +end + + +function population_factor(Ep::Vector{Float64}, T::Real)::Vector{Float64} + np = similar(Ep) + if isapprox(T,0.0,atol=PREC) + np[1] = 1.0 + np[2:end] .= 0.0 + return np + else + Z = partition_function(Ep, T) + @inbounds for i in eachindex(np) + np[i] = exp( -Ep[i]/(kB*T) ) / Z + end + return np + end +end + + +function thermal_average(; Ep::Vector{Float64}, Vp::Matrix{ComplexF64}, + op::Matrix{ComplexF64}, T::Real, mode::Function=real)::Float64 + tav::ComplexF64 = 0.0 + np = population_factor(Ep, T) + @views @inbounds for i in eachindex(Ep) + if isapprox(np[i], 0.0, atol=PREC) + continue + else + tav += dot(Vp[:, i], op, Vp[:, i]) * np[i] + end + end + return mode(tav) +end \ No newline at end of file