From 89d54982c43cd153ec0588fe4a258e9f7620f722 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 11:43:16 +0100 Subject: [PATCH 01/12] initial minimal commit single ion and cef matrix implementations --- .gitignore | 3 +- Manifest.toml | 303 ++++++++++++++++++++++++++++++++++++++++++ Project.toml | 18 ++- README.md | 4 +- src/CEF.jl | 47 ++++++- src/cef_system.jl | 327 ++++++++++++++++++++++++++++++++++++++++++++++ src/cef_utils.jl | 64 +++++++++ src/single_ion.jl | 217 ++++++++++++++++++++++++++++++ 8 files changed, 976 insertions(+), 7 deletions(-) create mode 100644 Manifest.toml create mode 100644 src/cef_system.jl create mode 100644 src/cef_utils.jl create mode 100644 src/single_ion.jl 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/src/CEF.jl b/src/CEF.jl index fbb0c88..27ddafc 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -1,5 +1,48 @@ 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_system +export cef_diagonalization +export cef_hamiltonian +export stevens_EO, stevens_O + +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..36e2c35 --- /dev/null +++ b/src/cef_system.jl @@ -0,0 +1,327 @@ +Base.@kwdef mutable struct cef_system + ion::mag_ion + cefparams::DataFrame + B::Vector{Float64} + H::HERMITIANC64 + E::Vector{Float64} + V::Matrix{ComplexF64} +end + + +function Base.show(io::IO, ::MIME"text/plain", cefsys::cef_system) + printstyled(io, "CEF matrix diagonalization results.\n\n", color=:underline, bold=true) + display(cefsys.ion) + println() + println("CEF parameters in (meV):") + println(cefsys.cefparams) + println() + println(io,"Diagonal g-tensor [gxx, gyy, gzz]: $(cefsys.ion.g).\n") + println(io,"External magnetic field in (Tesla) [Bx, By, Bz]: $(cefsys.B)\n") + println(io,"CEF energy levels in (meV) and in (Kelvin):") + E=cefsys.E + for i in eachindex(E) + Emev = round(E[i], digits=SDIG) + EK = round(E[i]/meV_per_K, digits=SDIG) + Es = Emev, EK + println(io,join(Es, ",\t")) + end + return nothing +end + + +@doc raw""" + cef_diagonalization(ion::mag_ion, bfactors::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::cef_system + +Display the CEF matrix diagonalization results, the effect of an external magnetic field `B` is modeled with a `Jz` (`Sz`) operator. +Returns `cef_system` struct with the result of the calculation. +Two calculation `methods` are implemented. `:EO`, which uses the efficient algorithm due to Rudowicz and Ryabov. +And `:O`, which calculates the Stevens operators via a lookup table. +""" +function cef_diagonalization(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::cef_system + cef_matrix = cef_hamiltonian(ion, cefparams; B=B, method=method) + @assert ishermitian(cef_matrix) + E,V = eigen(cef_matrix) + E .-= minimum(E) + return cef_system( + ion,cefparams,B,cef_matrix,E,V + ) +end + + +@doc raw""" + cef_hamiltonian(ion::mag_ion, bfactors::DataFrame; B::Vector{<:Real}=zeros(Float64, 3))::HERMITIANC64 + +Compute full Hamiltonian matrix given CEF parameters. +The matrix is outputted in the basis of CEF |J, mJ> wavefunctions with +mJ=-J, -J+1, ... J-1, J. + +If a magnetic field `B` is included, a Zeeman term is added to the Hamiltonian +where the g-factor of the ion `ion.g` is taken into account. +""" +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 + + +@doc raw""" + cef_hamiltonian(ion::mag_ion, D::Real, E::Real; B::Vector{<:Real}=zeros(Float64, 3))::HERMITIANC64 + +Compute full Hamiltonian matrix given effective anisotropy parameters `D` and `E`. + +The Hamiltonian is given by: D Jz^2 + E/2(Jp^2 + Jm^2). +The matrix is outputted in the basis of CEF |J, mJ> wavefunctions with +mJ=-J, -J+1, ... J-1, J. + +If a magnetic field `B` is included, a Zeeman term is added to the Hamiltonian +where the g-factor of the ion `ion.g` is taken into account. +""" +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 + + +@doc raw""" + stevens_O(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 + +Tabulated version of the Stevens EO by Martin Rotter in the McPhase manual +These resulting matrices for l in [2, 4, 6] and m = -l:1:l for J in 0.5:0.5:7.5 +are equivalent to those generated programmatically with the algorithm of +Ryabov and Rudowicz +""" +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..9bdd4e4 --- /dev/null +++ b/src/cef_utils.jl @@ -0,0 +1,64 @@ +@doc raw""" + blm_dframe(blm_dict::Dict{String, <:Real})::DataFrame + +Given a dictionary of Stevens coefficients of the form Blm -> Value, return +a DataFrame with equivalent information. +""" +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..9cbd039 --- /dev/null +++ b/src/single_ion.jl @@ -0,0 +1,217 @@ +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 + ], +) + + +@doc raw""" + single_ion(ion::String) + +Given a magnetic ion name, generate a `mag_ion` struct that contains +the following data: + +- Total angular momentum quantum number J, isotropic Landé-factor gJ +- Stevens geometric factors θ_l, α, β, γ +- Expectation value of radial wave functions , l={2, 4, 6} +- Magnetic form factor coefficients + +Currently only trivalent rare-earth ions are supported. +See `custom_spin` if other ions are required. + +Coefficients taken from Rotter Bauer (2010), +the McPhase online manual: https://mcphase.github.io/webpage/manual/node129.html +and Brown's coefficients https://www.ill.eu/sites/ccsl/ffacts/ffactnode1.html +""" +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 + + +@doc raw""" + spin_operators(J::Float64, a::String)::Matrix{ComplexF64} + +Explicit matrices for spin operators of given total angular momentum +quantum number `J`. +The string `a` is one of either `x`, `y`, `z` for the Cartesian spin-operators. +If `a` is either `+` or `-` explicit matrices for the ladder operators +are returned. +""" +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 From 9e65cccb44f1810d28094912eb72f957ec26237a Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 11:57:48 +0100 Subject: [PATCH 02/12] statistical quantities (thermal) --- src/CEF.jl | 5 ++++ src/thermodynamical_quantities.jl | 42 +++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+) create mode 100644 src/thermodynamical_quantities.jl diff --git a/src/CEF.jl b/src/CEF.jl index 27ddafc..69d0a32 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -45,4 +45,9 @@ export cef_diagonalization export cef_hamiltonian export stevens_EO, stevens_O +include("./thermodynamical_quantities.jl") +export partition_function +export population_factor +export thermal_average + 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 From 4fe5d4e83430b7b668a571e06a6008c4d0f33533 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 12:27:06 +0100 Subject: [PATCH 03/12] entropy and heat capacity calculation --- src/CEF.jl | 4 +++ src/cef_entropy.jl | 84 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+) create mode 100644 src/cef_entropy.jl diff --git a/src/CEF.jl b/src/CEF.jl index 69d0a32..8760a64 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -50,4 +50,8 @@ export partition_function export population_factor export thermal_average +include("./cef_entropy.jl") +export cef_entropy! +export cef_entropy_speclevels! + 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..f07816c --- /dev/null +++ b/src/cef_entropy.jl @@ -0,0 +1,84 @@ +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 + + +@doc raw""" + cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::String="SI")::Nothing + +Calculate the temperature-dependence of the magnetic entropy and specific +heat capacity of a magnetic system consisting of non-interacting magnetic ions. +The details of the CEF system are given in `cefsys`, which is a `cef_system` struct. + +The calculation is performed on a `DataFrame` `dfcalc`. +`dfcalc` must have columns `[:T, :Bx, :By, :Bz]`. + +The entropy and heat capacity are calculated in `:SI` units (J/mol/K). +""" +function cef_entropy!(cefsys::cef_system, 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(cefsys.ion,cefsys.cefparams;B=B,method=method)) + E .-= minimum(E) + :HC_CALC=mag_heatcap(E,:T)*convfac + end + dfcalc[:,:SM_CALC] = mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]) + return nothing +end + + +@doc raw""" + cef_entropy_speclevels!(single_ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; levels::UnitRange=1:4, units::String="SI")::Nothing + +Calculate the temperature-dependence of the magnetic entropy and specific +heat capacity of a magnetic system consisting of magnetic ions of type `single_ion`. +The CEF parameters are given in the `cefparams` `DataFrame`. +The calculation is performed on the `DataFrame` `dfcalc`. + +`dfcalc` must have column `[:T, :Bx, :By, :Bz]`. + +`levels` is a step range that specifies the indices of the contributing crystal +field energy levels. + +As per the default it is `1:4` which means that the first +4 levels contribute to the entropy of the system. + +The entropy and heat capacity are calculated in `SI` units [J/mol/K]. +""" +function cef_entropy_speclevels!(cefsys::cef_system, 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(cefsys.ion,cefsys.cefparams;B=B,method=method)) + E .-= minimum(E) + E=E[levels] + :HC_CALC=mag_heatcap(E,:T)*convfac + end + dfcalc[:,:SM_CALC] = mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]) + return nothing +end \ No newline at end of file From 61727bb8b49a9f4325efcddf87133a86844e9825 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 12:29:05 +0100 Subject: [PATCH 04/12] clean docs --- src/cef_entropy.jl | 32 +------------------------------- src/cef_system.jl | 46 ---------------------------------------------- src/cef_utils.jl | 6 ------ src/single_ion.jl | 27 --------------------------- 4 files changed, 1 insertion(+), 110 deletions(-) diff --git a/src/cef_entropy.jl b/src/cef_entropy.jl index f07816c..04a2cb4 100644 --- a/src/cef_entropy.jl +++ b/src/cef_entropy.jl @@ -2,7 +2,7 @@ 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]" + @error "Units not supported in heat capacity calculations. Use SI, (J/K/mol)" end end @@ -23,18 +23,6 @@ function mag_entropy(HC::Vector{Float64}, T::Vector{Float64})::Vector{Float64} end -@doc raw""" - cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::String="SI")::Nothing - -Calculate the temperature-dependence of the magnetic entropy and specific -heat capacity of a magnetic system consisting of non-interacting magnetic ions. -The details of the CEF system are given in `cefsys`, which is a `cef_system` struct. - -The calculation is performed on a `DataFrame` `dfcalc`. -`dfcalc` must have columns `[:T, :Bx, :By, :Bz]`. - -The entropy and heat capacity are calculated in `:SI` units (J/mol/K). -""" function cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:SI, method::Symbol=:EO)::Nothing convfac = hc_units(units) @eachrow! dfcalc begin @@ -49,24 +37,6 @@ function cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:SI, end -@doc raw""" - cef_entropy_speclevels!(single_ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; levels::UnitRange=1:4, units::String="SI")::Nothing - -Calculate the temperature-dependence of the magnetic entropy and specific -heat capacity of a magnetic system consisting of magnetic ions of type `single_ion`. -The CEF parameters are given in the `cefparams` `DataFrame`. -The calculation is performed on the `DataFrame` `dfcalc`. - -`dfcalc` must have column `[:T, :Bx, :By, :Bz]`. - -`levels` is a step range that specifies the indices of the contributing crystal -field energy levels. - -As per the default it is `1:4` which means that the first -4 levels contribute to the entropy of the system. - -The entropy and heat capacity are calculated in `SI` units [J/mol/K]. -""" function cef_entropy_speclevels!(cefsys::cef_system, dfcalc::DataFrame; levels::UnitRange=1:4, units::Symbol=:SI, method::Symbol=:EO)::Nothing # only levels specified contribute (2J+1 levels total) diff --git a/src/cef_system.jl b/src/cef_system.jl index 36e2c35..96b83e3 100644 --- a/src/cef_system.jl +++ b/src/cef_system.jl @@ -29,14 +29,6 @@ function Base.show(io::IO, ::MIME"text/plain", cefsys::cef_system) end -@doc raw""" - cef_diagonalization(ion::mag_ion, bfactors::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::cef_system - -Display the CEF matrix diagonalization results, the effect of an external magnetic field `B` is modeled with a `Jz` (`Sz`) operator. -Returns `cef_system` struct with the result of the calculation. -Two calculation `methods` are implemented. `:EO`, which uses the efficient algorithm due to Rudowicz and Ryabov. -And `:O`, which calculates the Stevens operators via a lookup table. -""" function cef_diagonalization(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::cef_system cef_matrix = cef_hamiltonian(ion, cefparams; B=B, method=method) @assert ishermitian(cef_matrix) @@ -48,16 +40,6 @@ function cef_diagonalization(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Rea end -@doc raw""" - cef_hamiltonian(ion::mag_ion, bfactors::DataFrame; B::Vector{<:Real}=zeros(Float64, 3))::HERMITIANC64 - -Compute full Hamiltonian matrix given CEF parameters. -The matrix is outputted in the basis of CEF |J, mJ> wavefunctions with -mJ=-J, -J+1, ... J-1, J. - -If a magnetic field `B` is included, a Zeeman term is added to the Hamiltonian -where the g-factor of the ion `ion.g` is taken into account. -""" 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) @@ -67,18 +49,6 @@ function cef_hamiltonian(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=z end -@doc raw""" - cef_hamiltonian(ion::mag_ion, D::Real, E::Real; B::Vector{<:Real}=zeros(Float64, 3))::HERMITIANC64 - -Compute full Hamiltonian matrix given effective anisotropy parameters `D` and `E`. - -The Hamiltonian is given by: D Jz^2 + E/2(Jp^2 + Jm^2). -The matrix is outputted in the basis of CEF |J, mJ> wavefunctions with -mJ=-J, -J+1, ... J-1, J. - -If a magnetic field `B` is included, a Zeeman term is added to the Hamiltonian -where the g-factor of the ion `ion.g` is taken into account. -""" 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)) @@ -161,14 +131,6 @@ function ryabov_clm(l::Int, m::Int)::Float64 end -@doc raw""" - stevens_EO(J::Real, l::Int, m::Int)::HERMITIANC64 - -Returns the explicit matrix for the `m` extended Stevens operator of rank `l` -for a given total angular momentum quantum number `J`. -A maximum rank of `l=7` is supported. -`m` is in the range `-l:1:l`. -""" function stevens_EO(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 T = ion.Jp^l # T^l_l for _ in l-1:-1:abs(m) # Eqn (1 and 2) of Ryabov (1999), see Stoll EasySpin @@ -185,14 +147,6 @@ function stevens_EO(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 end -@doc raw""" - stevens_O(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 - -Tabulated version of the Stevens EO by Martin Rotter in the McPhase manual -These resulting matrices for l in [2, 4, 6] and m = -l:1:l for J in 0.5:0.5:7.5 -are equivalent to those generated programmatically with the algorithm of -Ryabov and Rudowicz -""" function stevens_O(ion::mag_ion, l::Int, m::Int)::HERMITIANC64 J=ion.J m_dim = Int(2*J+1) diff --git a/src/cef_utils.jl b/src/cef_utils.jl index 9bdd4e4..3970fc3 100644 --- a/src/cef_utils.jl +++ b/src/cef_utils.jl @@ -1,9 +1,3 @@ -@doc raw""" - blm_dframe(blm_dict::Dict{String, <:Real})::DataFrame - -Given a dictionary of Stevens coefficients of the form Blm -> Value, return -a DataFrame with equivalent information. -""" function blm_dframe(blm_dict::Dict{String, <:Real})::DataFrame l = Int[] m = Int[] diff --git a/src/single_ion.jl b/src/single_ion.jl index 9cbd039..f5c2600 100644 --- a/src/single_ion.jl +++ b/src/single_ion.jl @@ -79,24 +79,6 @@ const re_hundsrules = Dict( ) -@doc raw""" - single_ion(ion::String) - -Given a magnetic ion name, generate a `mag_ion` struct that contains -the following data: - -- Total angular momentum quantum number J, isotropic Landé-factor gJ -- Stevens geometric factors θ_l, α, β, γ -- Expectation value of radial wave functions , l={2, 4, 6} -- Magnetic form factor coefficients - -Currently only trivalent rare-earth ions are supported. -See `custom_spin` if other ions are required. - -Coefficients taken from Rotter Bauer (2010), -the McPhase online manual: https://mcphase.github.io/webpage/manual/node129.html -and Brown's coefficients https://www.ill.eu/sites/ccsl/ffacts/ffactnode1.html -""" function single_ion(ion::String) try J, g, @@ -175,15 +157,6 @@ function Base.show(io::IO, ::MIME"text/plain", ion::mag_ion) end -@doc raw""" - spin_operators(J::Float64, a::String)::Matrix{ComplexF64} - -Explicit matrices for spin operators of given total angular momentum -quantum number `J`. -The string `a` is one of either `x`, `y`, `z` for the Cartesian spin-operators. -If `a` is either `+` or `-` explicit matrices for the ladder operators -are returned. -""" function spin_operators(J::Float64, a::String)::Matrix{ComplexF64} mJ = -J:1:J if isequal(a, "x") From bf801810e2d67f51bf3c54bbc601afcc45167fa1 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 14:26:09 +0100 Subject: [PATCH 05/12] add magnetization calculations --- Manifest.toml | 123 ++++++++++++++++++++++++++++++++++++++- Project.toml | 2 + src/CEF.jl | 13 +++-- src/cef_magnetization.jl | 78 +++++++++++++++++++++++++ 4 files changed, 211 insertions(+), 5 deletions(-) create mode 100644 src/cef_magnetization.jl diff --git a/Manifest.toml b/Manifest.toml index 14d21c1..fa49f6d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,11 @@ julia_version = "1.11.1" manifest_format = "2.0" -project_hash = "0978e35ac3c4ed2f2567f101098e29d371b1f36c" +project_hash = "5e2bd8eb2b56513d1bea211b45c3088377dd09e5" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.2" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -70,6 +74,15 @@ deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" version = "1.11.0" +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" +version = "1.11.0" + [[deps.Future]] deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" @@ -88,6 +101,12 @@ version = "1.4.3" ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +[[deps.IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] +git-tree-sha1 = "0f14a5456bdc6b9731a5682f439a672750a09e48" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2025.0.4+0" + [[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -103,11 +122,47 @@ git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" uuid = "82899510-4779-5014-852e-03e436cf321d" version = "1.0.0" +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.7.0" + [[deps.LaTeXStrings]] git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" version = "1.4.0" +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" +version = "1.11.0" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.4" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "8.6.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +version = "1.11.0" + +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.7.2+0" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.11.0+1" + [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" version = "1.11.0" @@ -117,6 +172,22 @@ deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" version = "1.11.0" +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +version = "1.11.0" + +[[deps.MKL]] +deps = ["Artifacts", "Libdl", "LinearAlgebra", "Logging", "MKL_jll"] +git-tree-sha1 = "39fd578c295085e77aedec501149a63f9ab0c0b6" +uuid = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +version = "0.8.0" + +[[deps.MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] +git-tree-sha1 = "5de60bc6cb3899cd318d80d627560fae2e2d99ae" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2025.0.1+1" + [[deps.MacroTools]] git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" @@ -127,12 +198,25 @@ deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" version = "1.11.0" +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.6+0" + [[deps.Missings]] deps = ["DataAPI"] git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" version = "1.2.0" +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2023.12.12" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + [[deps.OffsetArrays]] git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -154,6 +238,17 @@ git-tree-sha1 = "cc4054e898b852042d7b503313f7ad03de99c3dd" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.8.0" +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.11.0" + + [deps.Pkg.extensions] + REPLExt = "REPL" + + [deps.Pkg.weakdeps] + REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + [[deps.PooledArrays]] deps = ["DataAPI", "Future"] git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3" @@ -269,6 +364,11 @@ git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" version = "1.12.0" +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + [[deps.Trapz]] git-tree-sha1 = "79eb0ed763084a3e7de81fe1838379ac6a23b6a0" uuid = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" @@ -297,7 +397,28 @@ version = "1.22.0" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+1" + [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" version = "5.11.0+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.59.0+0" + +[[deps.oneTBB_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "d5a767a3bb77135a99e433afe0eb14cd7f6914c3" +uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" +version = "2022.0.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+2" diff --git a/Project.toml b/Project.toml index 4368344..5c33853 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.0-DEV" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -16,6 +17,7 @@ Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" DataFrames = "1.7.0" DataFramesMeta = "0.15.4" LinearAlgebra = "1.11.0" +MKL = "0.8.0" OffsetArrays = "1.15.0" StaticArrays = "1.9.13" Statistics = "1.11.1" diff --git a/src/CEF.jl b/src/CEF.jl index 8760a64..09e99de 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -7,6 +7,7 @@ using OffsetArrays using StaticArrays using Statistics using Trapz +using MKL const PREC::Float64 = 1.0e-7 # for degeneracy calculations const SDIG::Int64 = 7 # for numerical cutoffs @@ -28,6 +29,10 @@ 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("./cef_utils.jl") +export blm_dframe +export get_alm! + include("./single_ion.jl") export single_ion export custom_ion @@ -35,10 +40,6 @@ 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_system export cef_diagonalization @@ -54,4 +55,8 @@ include("./cef_entropy.jl") export cef_entropy! export cef_entropy_speclevels! +include("./cef_magnetization.jl") +export cef_magneticmoment_crystal! +export cef_magneticmoment_powder! + 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..52ae564 --- /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!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) + unit_factor = mag_units(units) + spinops = [cefsys.ion.Jx, cefsys.ion.Jy, cefsys.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(cefsys.ion,cefsys.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 = norm(cefsys.ion.g .* spin_expval)*unit_factor + :Mx_CALC = thermal_average(Ep=E,Vp=V,op=spinops[1],T=:T,mode=mode)*cefsys.ion.g[1]*unit_factor + :My_CALC = thermal_average(Ep=E,Vp=V,op=spinops[2],T=:T,mode=mode)*cefsys.ion.g[2]*unit_factor + :Mz_CALC = thermal_average(Ep=E,Vp=V,op=spinops[3],T=:T,mode=mode)*cefsys.ion.g[3]*unit_factor + end + return nothing +end + + +function cef_magneticmoment_powder!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) + unit_factor = mag_units(units) + spinops = [cefsys.ion.Jx, cefsys.ion.Jy, cefsys.ion.Jz] + @eachrow! dfcalc begin + @newcol :M_CALC::Vector{Float64} + E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[:B,0.0,0.0],method=method)) + E .-= minimum(E) + MX = calc_magmom(cefsys.ion.g,spinops .* [1.0, 0.0, 0.0],E,V,:T,mode) + + E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[0.0,:B,0.0],method=method)) + E .-= minimum(E) + MY = calc_magmom(cefsys.ion.g,spinops .* [0.0, 1.0, 0.0],E,V,:T,mode) + + E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[0.0,0.0,:B],method=method)) + E .-= minimum(E) + MZ = calc_magmom(cefsys.ion.g,spinops .* [0.0, 0.0, 1.0],E,V,:T,mode) + + :M_CALC = ((MX + MY + MZ) / 3.0) * unit_factor + end + return nothing +end \ No newline at end of file From 43448fcccd79fc24c2d2d5ab2b3304bff97252ab Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 19:28:52 +0100 Subject: [PATCH 06/12] rounding to significant figures (7) --- Manifest.toml | 123 +-------------------------------------- Project.toml | 2 - src/CEF.jl | 9 ++- src/cef_entropy.jl | 8 +-- src/cef_magnetization.jl | 10 ++-- 5 files changed, 14 insertions(+), 138 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index fa49f6d..14d21c1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,11 +2,7 @@ julia_version = "1.11.1" manifest_format = "2.0" -project_hash = "5e2bd8eb2b56513d1bea211b45c3088377dd09e5" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.2" +project_hash = "0978e35ac3c4ed2f2567f101098e29d371b1f36c" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -74,15 +70,6 @@ deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" version = "1.11.0" -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" -version = "1.11.0" - [[deps.Future]] deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" @@ -101,12 +88,6 @@ version = "1.4.3" ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -[[deps.IntelOpenMP_jll]] -deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] -git-tree-sha1 = "0f14a5456bdc6b9731a5682f439a672750a09e48" -uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2025.0.4+0" - [[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -122,47 +103,11 @@ git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" uuid = "82899510-4779-5014-852e-03e436cf321d" version = "1.0.0" -[[deps.JLLWrappers]] -deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.7.0" - [[deps.LaTeXStrings]] git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" version = "1.4.0" -[[deps.LazyArtifacts]] -deps = ["Artifacts", "Pkg"] -uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" -version = "1.11.0" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.6.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -version = "1.11.0" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.7.2+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" - [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" version = "1.11.0" @@ -172,22 +117,6 @@ deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" version = "1.11.0" -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -version = "1.11.0" - -[[deps.MKL]] -deps = ["Artifacts", "Libdl", "LinearAlgebra", "Logging", "MKL_jll"] -git-tree-sha1 = "39fd578c295085e77aedec501149a63f9ab0c0b6" -uuid = "33e6dc65-8f57-5167-99aa-e5a354878fb2" -version = "0.8.0" - -[[deps.MKL_jll]] -deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] -git-tree-sha1 = "5de60bc6cb3899cd318d80d627560fae2e2d99ae" -uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2025.0.1+1" - [[deps.MacroTools]] git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" @@ -198,25 +127,12 @@ deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" version = "1.11.0" -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.6+0" - [[deps.Missings]] deps = ["DataAPI"] git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" version = "1.2.0" -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.12.12" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - [[deps.OffsetArrays]] git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -238,17 +154,6 @@ git-tree-sha1 = "cc4054e898b852042d7b503313f7ad03de99c3dd" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.8.0" -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.11.0" - - [deps.Pkg.extensions] - REPLExt = "REPL" - - [deps.Pkg.weakdeps] - REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - [[deps.PooledArrays]] deps = ["DataAPI", "Future"] git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3" @@ -364,11 +269,6 @@ git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" version = "1.12.0" -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - [[deps.Trapz]] git-tree-sha1 = "79eb0ed763084a3e7de81fe1838379ac6a23b6a0" uuid = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" @@ -397,28 +297,7 @@ version = "1.22.0" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+1" - [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" version = "5.11.0+0" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.59.0+0" - -[[deps.oneTBB_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "d5a767a3bb77135a99e433afe0eb14cd7f6914c3" -uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" -version = "2022.0.0+0" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+2" diff --git a/Project.toml b/Project.toml index 5c33853..4368344 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,6 @@ version = "0.1.0-DEV" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -17,7 +16,6 @@ Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" DataFrames = "1.7.0" DataFramesMeta = "0.15.4" LinearAlgebra = "1.11.0" -MKL = "0.8.0" OffsetArrays = "1.15.0" StaticArrays = "1.9.13" Statistics = "1.11.1" diff --git a/src/CEF.jl b/src/CEF.jl index 09e99de..e00cf53 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -7,7 +7,6 @@ using OffsetArrays using StaticArrays using Statistics using Trapz -using MKL const PREC::Float64 = 1.0e-7 # for degeneracy calculations const SDIG::Int64 = 7 # for numerical cutoffs @@ -29,10 +28,6 @@ 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("./cef_utils.jl") -export blm_dframe -export get_alm! - include("./single_ion.jl") export single_ion export custom_ion @@ -40,6 +35,10 @@ 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_system export cef_diagonalization diff --git a/src/cef_entropy.jl b/src/cef_entropy.jl index 04a2cb4..16f93e7 100644 --- a/src/cef_entropy.jl +++ b/src/cef_entropy.jl @@ -30,9 +30,9 @@ function cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:SI, B=[:Bx,:By,:Bz] E=eigvals(cef_hamiltonian(cefsys.ion,cefsys.cefparams;B=B,method=method)) E .-= minimum(E) - :HC_CALC=mag_heatcap(E,:T)*convfac + :HC_CALC=round(mag_heatcap(E,:T)*convfac,digits=SDIG) end - dfcalc[:,:SM_CALC] = mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]) + dfcalc[:,:SM_CALC]=round(mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]),digits=SDIG) return nothing end @@ -47,8 +47,8 @@ function cef_entropy_speclevels!(cefsys::cef_system, dfcalc::DataFrame; E=eigvals(cef_hamiltonian(cefsys.ion,cefsys.cefparams;B=B,method=method)) E .-= minimum(E) E=E[levels] - :HC_CALC=mag_heatcap(E,:T)*convfac + :HC_CALC=round(mag_heatcap(E,:T)*convfac,digits=SDIG) end - dfcalc[:,:SM_CALC] = mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]) + 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 index 52ae564..9403770 100644 --- a/src/cef_magnetization.jl +++ b/src/cef_magnetization.jl @@ -46,10 +46,10 @@ function cef_magneticmoment_crystal!(cefsys::cef_system, dfcalc::DataFrame; unit end spin_expval[i] = thermal_average(Ep=E,Vp=V,op=spinops[i],T=:T,mode=mode) end - :Mp_CALC = norm(cefsys.ion.g .* spin_expval)*unit_factor - :Mx_CALC = thermal_average(Ep=E,Vp=V,op=spinops[1],T=:T,mode=mode)*cefsys.ion.g[1]*unit_factor - :My_CALC = thermal_average(Ep=E,Vp=V,op=spinops[2],T=:T,mode=mode)*cefsys.ion.g[2]*unit_factor - :Mz_CALC = thermal_average(Ep=E,Vp=V,op=spinops[3],T=:T,mode=mode)*cefsys.ion.g[3]*unit_factor + :Mp_CALC=round(norm(cefsys.ion.g .* spin_expval)*unit_factor,digits=SDIG) + :Mx_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[1],T=:T,mode=mode)*cefsys.ion.g[1]*unit_factor,digits=SDIG) + :My_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[2],T=:T,mode=mode)*cefsys.ion.g[2]*unit_factor,digits=SDIG) + :Mz_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[3],T=:T,mode=mode)*cefsys.ion.g[3]*unit_factor,digits=SDIG) end return nothing end @@ -72,7 +72,7 @@ function cef_magneticmoment_powder!(cefsys::cef_system, dfcalc::DataFrame; units E .-= minimum(E) MZ = calc_magmom(cefsys.ion.g,spinops .* [0.0, 0.0, 1.0],E,V,:T,mode) - :M_CALC = ((MX + MY + MZ) / 3.0) * unit_factor + :M_CALC=round(((MX + MY + MZ) / 3.0) * unit_factor,digits=SDIG) end return nothing end \ No newline at end of file From 468d74ed6ad8800ed6bb2a2508afb75213939e5c Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 19:49:40 +0100 Subject: [PATCH 07/12] rollback implementation of cefsys struct --- src/CEF.jl | 5 ++++- src/cef_entropy.jl | 9 ++++---- src/cef_magnetization.jl | 30 ++++++++++++------------- src/cef_system.jl | 48 ++++++++++++++-------------------------- 4 files changed, 40 insertions(+), 52 deletions(-) diff --git a/src/CEF.jl b/src/CEF.jl index e00cf53..15bc17f 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -40,7 +40,6 @@ export blm_dframe export get_alm! include("./cef_system.jl") -export cef_system export cef_diagonalization export cef_hamiltonian export stevens_EO, stevens_O @@ -58,4 +57,8 @@ include("./cef_magnetization.jl") export cef_magneticmoment_crystal! export cef_magneticmoment_powder! +include("./cef_susceptibility.jl") +export cef_susceptibility_crystal! +export cef_susceptibility_powder! + end \ No newline at end of file diff --git a/src/cef_entropy.jl b/src/cef_entropy.jl index 16f93e7..82995d4 100644 --- a/src/cef_entropy.jl +++ b/src/cef_entropy.jl @@ -23,12 +23,12 @@ function mag_entropy(HC::Vector{Float64}, T::Vector{Float64})::Vector{Float64} end -function cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:SI, method::Symbol=:EO)::Nothing +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(cefsys.ion,cefsys.cefparams;B=B,method=method)) + E=eigvals(cef_hamiltonian(ion,cefparams;B=B,method=method)) E .-= minimum(E) :HC_CALC=round(mag_heatcap(E,:T)*convfac,digits=SDIG) end @@ -37,14 +37,13 @@ function cef_entropy!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:SI, end -function cef_entropy_speclevels!(cefsys::cef_system, dfcalc::DataFrame; - levels::UnitRange=1:4, units::Symbol=:SI, method::Symbol=:EO)::Nothing +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(cefsys.ion,cefsys.cefparams;B=B,method=method)) + 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) diff --git a/src/cef_magnetization.jl b/src/cef_magnetization.jl index 9403770..45558a5 100644 --- a/src/cef_magnetization.jl +++ b/src/cef_magnetization.jl @@ -23,16 +23,16 @@ function calc_magmom(g::MVEC{3}, spinops::Vector{Matrix{ComplexF64}}, Ep::Vector end -function cef_magneticmoment_crystal!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) +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 = [cefsys.ion.Jx, cefsys.ion.Jy, cefsys.ion.Jz] + 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(cefsys.ion,cefsys.cefparams;B=extfield,method=method)) + E, V = eigen(cef_hamiltonian(ion,cefparams;B=extfield,method=method)) E .-= minimum(E) if iszero(extfield) spin_proj = spinops @@ -46,31 +46,31 @@ function cef_magneticmoment_crystal!(cefsys::cef_system, dfcalc::DataFrame; unit end spin_expval[i] = thermal_average(Ep=E,Vp=V,op=spinops[i],T=:T,mode=mode) end - :Mp_CALC=round(norm(cefsys.ion.g .* spin_expval)*unit_factor,digits=SDIG) - :Mx_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[1],T=:T,mode=mode)*cefsys.ion.g[1]*unit_factor,digits=SDIG) - :My_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[2],T=:T,mode=mode)*cefsys.ion.g[2]*unit_factor,digits=SDIG) - :Mz_CALC=round(thermal_average(Ep=E,Vp=V,op=spinops[3],T=:T,mode=mode)*cefsys.ion.g[3]*unit_factor,digits=SDIG) + :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!(cefsys::cef_system, dfcalc::DataFrame; units::Symbol=:ATOMIC, method::Symbol=:EO, mode::Function=real) +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 = [cefsys.ion.Jx, cefsys.ion.Jy, cefsys.ion.Jz] + spinops = [ion.Jx,ion.Jy,ion.Jz] @eachrow! dfcalc begin @newcol :M_CALC::Vector{Float64} - E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[:B,0.0,0.0],method=method)) + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[:B,0.0,0.0],method=method)) E .-= minimum(E) - MX = calc_magmom(cefsys.ion.g,spinops .* [1.0, 0.0, 0.0],E,V,:T,mode) + MX = calc_magmom(ion.g,spinops .* [1.0, 0.0, 0.0],E,V,:T,mode) - E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[0.0,:B,0.0],method=method)) + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[0.0,:B,0.0],method=method)) E .-= minimum(E) - MY = calc_magmom(cefsys.ion.g,spinops .* [0.0, 1.0, 0.0],E,V,:T,mode) + MY = calc_magmom(ion.g,spinops .* [0.0, 1.0, 0.0],E,V,:T,mode) - E, V = eigen(cef_hamiltonian(cefsys.ion,cefsys.cefparams; B=[0.0,0.0,:B],method=method)) + E, V = eigen(cef_hamiltonian(ion,cefparams; B=[0.0,0.0,:B],method=method)) E .-= minimum(E) - MZ = calc_magmom(cefsys.ion.g,spinops .* [0.0, 0.0, 1.0],E,V,:T,mode) + 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 diff --git a/src/cef_system.jl b/src/cef_system.jl index 96b83e3..3712f5a 100644 --- a/src/cef_system.jl +++ b/src/cef_system.jl @@ -1,45 +1,31 @@ -Base.@kwdef mutable struct cef_system - ion::mag_ion - cefparams::DataFrame - B::Vector{Float64} - H::HERMITIANC64 - E::Vector{Float64} - V::Matrix{ComplexF64} -end +@doc raw""" + cef_eigensystem(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3))::Nothing +Display CEF matrix diagonalization results. Useful in cases where a quick +view of the CEF energy spectrum is desired. -function Base.show(io::IO, ::MIME"text/plain", cefsys::cef_system) - printstyled(io, "CEF matrix diagonalization results.\n\n", color=:underline, bold=true) - display(cefsys.ion) - println() - println("CEF parameters in (meV):") - println(cefsys.cefparams) - println() - println(io,"Diagonal g-tensor [gxx, gyy, gzz]: $(cefsys.ion.g).\n") - println(io,"External magnetic field in (Tesla) [Bx, By, Bz]: $(cefsys.B)\n") - println(io,"CEF energy levels in (meV) and in (Kelvin):") - E=cefsys.E +Details of the magnetic ion, the CEF parameters and applied field are displayed. +""" +function cef_eigensystem(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 Kelvin:") for i in eachindex(E) Emev = round(E[i], digits=SDIG) EK = round(E[i]/meV_per_K, digits=SDIG) Es = Emev, EK - println(io,join(Es, ",\t")) + println(join(Es, ",\t")) end return nothing end -function cef_diagonalization(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::cef_system - cef_matrix = cef_hamiltonian(ion, cefparams; B=B, method=method) - @assert ishermitian(cef_matrix) - E,V = eigen(cef_matrix) - E .-= minimum(E) - return cef_system( - ion,cefparams,B,cef_matrix,E,V - ) -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) From 9989369e20dc07b134971c9cff5c38bd57463041 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 19:49:53 +0100 Subject: [PATCH 08/12] add susceptibility calculations --- src/cef_susceptibility.jl | 91 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 src/cef_susceptibility.jl diff --git a/src/cef_susceptibility.jl b/src/cef_susceptibility.jl new file mode 100644 index 0000000..33446b1 --- /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,bfactors;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 From 582d22ddd7c7e8507e350f10d3fccfe9a678702d Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 20:21:04 +0100 Subject: [PATCH 09/12] add neutron x section --- src/CEF.jl | 4 ++ src/cef_neutronxsection.jl | 144 +++++++++++++++++++++++++++++++++++++ src/cef_system.jl | 14 +--- 3 files changed, 151 insertions(+), 11 deletions(-) create mode 100644 src/cef_neutronxsection.jl diff --git a/src/CEF.jl b/src/CEF.jl index 15bc17f..2d6cb2c 100644 --- a/src/CEF.jl +++ b/src/CEF.jl @@ -61,4 +61,8 @@ 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_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_system.jl b/src/cef_system.jl index 3712f5a..3a73ac8 100644 --- a/src/cef_system.jl +++ b/src/cef_system.jl @@ -1,12 +1,4 @@ -@doc raw""" - cef_eigensystem(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3))::Nothing - -Display CEF matrix diagonalization results. Useful in cases where a quick -view of the CEF energy spectrum is desired. - -Details of the magnetic ion, the CEF parameters and applied field are displayed. -""" -function cef_eigensystem(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=zeros(Float64, 3), method::Symbol=:EO)::Nothing +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) @@ -14,8 +6,8 @@ function cef_eigensystem(ion::mag_ion, cefparams::DataFrame; B::Vector{<:Real}=z 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 Kelvin:") + 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) From 1c7cd65f107dd008ccb9deac2d3ce29e2140755c Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 20:30:43 +0100 Subject: [PATCH 10/12] add examples --- examples/hex_yb.jl | 68 ++++++++++++++++++++++++++++++++ examples/re_neutronformfactor.jl | 38 ++++++++++++++++++ src/cef_susceptibility.jl | 2 +- 3 files changed, 107 insertions(+), 1 deletion(-) create mode 100644 examples/hex_yb.jl create mode 100644 examples/re_neutronformfactor.jl diff --git a/examples/hex_yb.jl b/examples/hex_yb.jl new file mode 100644 index 0000000..6072e2c --- /dev/null +++ b/examples/hex_yb.jl @@ -0,0 +1,68 @@ +""" +Reproduction of results from the paper by Rotter and Bauer +""" + + +using CEF +using DataFrames +using StatsPlots + + +ion = single_ion("Yb3+") +bfactors = blm_dframe(Dict("B20"=>0.5622, "B40"=>1.6087e-5, + "B60"=>6.412e-7, "B66"=>-8.324e-6)) + + +function main() + + """ 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") + # chi_powd = 1/3*calc_grid_z.CHI_CALC + 2/3*calc_grid_x.CHI_CALC + # plot!(chi_plot, calc_grid_z.T, 1 ./ chi_powd, 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..6153f22 --- /dev/null +++ b/examples/re_neutronformfactor.jl @@ -0,0 +1,38 @@ +""" +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->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_susceptibility.jl b/src/cef_susceptibility.jl index 33446b1..e4f186c 100644 --- a/src/cef_susceptibility.jl +++ b/src/cef_susceptibility.jl @@ -78,7 +78,7 @@ 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,bfactors;method=method)) + E, V = eigen(cef_hamiltonian(ion,cefparams;method=method)) E .-= minimum(E) @eachrow! dfcalc begin @newcol :CHI_CALC::Vector{Float64} From 42bb93c6d54c7435057bab8aaa4808615b1b396a Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 20:32:44 +0100 Subject: [PATCH 11/12] fix broadcast round --- src/cef_entropy.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cef_entropy.jl b/src/cef_entropy.jl index 82995d4..a26e047 100644 --- a/src/cef_entropy.jl +++ b/src/cef_entropy.jl @@ -32,7 +32,7 @@ function cef_entropy!(ion::mag_ion, cefparams::DataFrame, dfcalc::DataFrame; uni 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) + dfcalc[:,:SM_CALC]=round.(mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]),digits=SDIG) return nothing end @@ -48,6 +48,6 @@ function cef_entropy_speclevels!(ion::mag_ion, cefparams::DataFrame, dfcalc::Dat 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) + dfcalc[:,:SM_CALC]=round.(mag_entropy(dfcalc[:,:HC_CALC],dfcalc[:,:T]),digits=SDIG) return nothing end \ No newline at end of file From a1d3368a379a2371a53b6996ef98370f040a6b78 Mon Sep 17 00:00:00 2001 From: abehersan Date: Mon, 3 Mar 2025 21:11:41 +0100 Subject: [PATCH 12/12] update examples --- examples/hex_yb.jl | 14 ++++---------- examples/re_neutronformfactor.jl | 3 ++- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/examples/hex_yb.jl b/examples/hex_yb.jl index 6072e2c..14408db 100644 --- a/examples/hex_yb.jl +++ b/examples/hex_yb.jl @@ -8,13 +8,12 @@ using DataFrames using StatsPlots -ion = single_ion("Yb3+") -bfactors = blm_dframe(Dict("B20"=>0.5622, "B40"=>1.6087e-5, - "B60"=>6.412e-7, "B66"=>-8.324e-6)) - - 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]") @@ -24,7 +23,6 @@ function main() @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 @@ -36,7 +34,6 @@ function main() @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) @@ -48,9 +45,6 @@ function main() 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") - # chi_powd = 1/3*calc_grid_z.CHI_CALC + 2/3*calc_grid_x.CHI_CALC - # plot!(chi_plot, calc_grid_z.T, 1 ./ chi_powd, label="Powder") - """ Specific heat capacity and magnetic entropy """ s_plot = plot(xlabel="Temperature [K]", ylabel="S, HC (J/mol/K)") diff --git a/examples/re_neutronformfactor.jl b/examples/re_neutronformfactor.jl index 6153f22..63967b6 100644 --- a/examples/re_neutronformfactor.jl +++ b/examples/re_neutronformfactor.jl @@ -11,6 +11,7 @@ using StatsPlots function main() + mag_ions = [ "Ce3+", "Pr3+", "Nd3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+" ] @@ -22,7 +23,7 @@ function main() for i in eachindex(mag_ions) ion = single_ion(mag_ions[i]) lbl = "Ion: $(ion.ion), J: $(ion.J)" - FSQUARED = abs2.(broadcast(q->dipolar_formfactor(ion, q), QS)) + 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)