Dev

.gitignore

@@ -1,2 +1,2 @@
 ./Manifest.toml
-./github
+./.github

src/CEF.jl

@@ -65,4 +65,7 @@ include("./cef_neutronxsection.jl")
 export cef_neutronxsection_crystal!
 export cef_neutronxsection_powder!
 
+include("./cef_rotation.jl")
+export rotate_blm
+
 end

src/cef_rotation.jl

@@ -0,0 +1,145 @@
+function rotation_unitary(J::Real, n::Vector{<:Real}, phi::Real)::Matrix{ComplexF64}
+    Jx = spin_operators(J, "x")
+    Jy = spin_operators(J, "y")
+    Jz = spin_operators(J, "z")
+    nnorm = normalize(n)
+    Jproj = sum(nnorm .* [Jx, Jy, Jz])
+    U = exp(-1im * phi * Jproj)
+    @assert isunitary(U)
+    return U
+end
+
+
+function ZYZ_rotmatrix(alpha::Real, beta::Real, gamma::Real)::Matrix{Float64}
+    a, b, g = map(Float64, [alpha, beta, gamma])
+    return Z_rot(a) * Y_rot(b) * Z_rot(g)
+end
+
+
+X_rot(theta::Float64)::Matrix{Float64} = [1 0 0;
+                                          0 cos(theta) -sin(theta);
+                                          0 sin(theta) cos(theta)]
+
+
+Y_rot(theta::Float64)::Matrix{Float64} = [cos(theta) 0 sin(theta);
+                                          0 1 0;
+                                          -sin(theta) 0 cos(theta)]
+
+
+Z_rot(theta::Float64)::Matrix{Float64} = [cos(theta) -sin(theta) 0;
+                                          sin(theta) cos(theta) 0;
+                                          0 0 1]
+
+
+function get_euler_angles(v::Vector{<:Real})::Tuple{Float64, Float64, Float64}
+    if isequal(zeros(Real, 3), v)
+        return (0.0, 0.0, 0.0)
+    end
+    v_norm = v / norm(v)
+    alpha = atan(v_norm[2], v_norm[1])# / pi * 180.0
+    beta = atan((v_norm[1]^2 + v_norm[2]^2), v_norm[3])# / pi * 180.0
+    gamma = 0.0
+    return (alpha, beta, gamma) # in radian
+end
+
+
+function wigner_D(l::Int, alpha::Real, beta::Real, gamma::Real)::Matrix{ComplexF64}
+    mdim = Int(2*l+1)
+    ml = -l:1:l
+    rotmat = Matrix{ComplexF64}(undef, (mdim, mdim))
+    for (idmp, mm) in enumerate(ml)
+        for (idm, mp) in enumerate(ml)
+            rotmat[idmp, idm] = (-1)^(mp - mm) * rotation_matrix_element(l, mm, mp, alpha, beta, gamma)
+        end
+    end
+    return round.(rotmat, digits=SDIG)
+end
+
+
+function rotation_matrix_element(l::Int, m::Int, mp::Int, alpha::Real, beta::Real, gamma::Real)::ComplexF64
+    return exp(-1.0im*mp*alpha) * small_d(l, mp, m, beta) * exp(-1.0im*m*gamma)
+end
+
+
+function small_d(l::Int, mp::Int, m::Int, beta::Real)::Float64
+    if iszero(beta)
+        return isequal(m, mp) * 1.0 # delta function d_m'm if beta=0
+    end
+    djmpm = 0.0
+    smin = Int(maximum([0, -(m+mp)]))
+    smax = Int(minimum([l-m, l-mp]))
+    for s in smin:1:smax
+        djmpm += (-1)^(l-m-s)*
+            ((cos(beta/2))^(2s+mp+m))*
+            ((sin(beta/2))^(2l-2s-m-mp))*
+            binomial(Int(l+m), Int(l-mp-s))*
+            binomial(Int(l-m), Int(s))
+    end
+    djmpm *= sqrt((factorial(Int(l+mp))*factorial(Int(l-mp)))/
+                  (factorial(Int(l+m))*factorial(Int(l-m))))
+    if isequal(djmpm, NaN) || isequal(djmpm, Inf)
+        err_message =
+        "Matrix element djmpm is NaN or Inf!\n"*
+        "Likely there's an issue in the inputted values of l, m or mp."
+        @error err_message
+    else
+        return djmpm
+    end
+end
+
+
+function rotate_blm(cefparams::DataFrame, alpha::Real, beta::Real, gamma::Real)::DataFrame
+    cefparams_rotated = DataFrame(B = Float64[], l = Int[], m = Int[])
+    ls = sort(collect(Set(cefparams[:, :l])))
+    for l in ls
+        S_matrix = rotate_stevens(l, alpha, beta, gamma)
+        cefparams_ori = zeros(Float64, Int(2l+1)) # original CEF parameters
+        cefparams_rot = zeros(ComplexF64, Int(2l+1)) # rotated CEF params (complex)
+        cefparams_res = zeros(Float64, Int(2l+1)) # rotated CEF parameters (real)
+        for (i, m) in enumerate(-l:1:l)
+            try
+                cefparams_ori[i] = cefparams[(cefparams.l .== l) .& (cefparams.m .== m), :B][1]
+            catch ex
+                if isa(ex, BoundsError)
+                    cefparams_ori[i] = 0.0
+                end
+            end
+        end
+        # S * cefparams where cefparams is a (2l+1) vector
+        cefparams_rot .= S_matrix * cefparams_ori
+        @assert norm(imag(cefparams_rot)) < PREC
+        cefparams_res .= real(cefparams_rot)
+        append!(cefparams_rotated, DataFrame("B"=>cefparams_res, "l"=>fill(l, Int(2l+1)), "m"=>-l:1:l))
+    end
+    return cefparams_rotated
+end
+
+
+function rotate_stevens(l::Int, alpha::Real, beta::Real, gamma::Real)::Matrix{ComplexF64}
+    return transpose(inv(Alm_matrix(l))) * wigner_D(l, alpha, beta, gamma)' * transpose(Alm_matrix(l))
+end
+
+
+function Alm_matrix(l::Int)::Matrix{ComplexF64}
+    alm_coeff = OffsetArray(SMatrix{6, 7}([
+        1           1/sqrt(2)       0               0               0               0           0;
+        sqrt(6)     1/2             1               0               0               0           0;
+        sqrt(10)    sqrt(10/3)      1/sqrt(3)       sqrt(2)         0               0           0;
+        2*sqrt(70)  sqrt(7/2)       sqrt(7)         1/sqrt(2)       2               2           0;
+        6*sqrt(14)  2*sqrt(21/5)    sqrt(3/5)       6*sqrt(2/5)     2/sqrt(5)       2*sqrt(2)   0;
+        4*sqrt(231) sqrt(22)        4*sqrt(11/5)    2*sqrt(11/5)    4*sqrt(11/6)    2/sqrt(3)   4;
+    ]), 1:6, 0:6)
+    a_matrix = OffsetArray(zeros(ComplexF64, (2l+1, 2l+1)), -l:l, -l:l)
+    for m in -l:1:l
+        if iszero(m)
+            a_matrix[m, m] = alm_coeff[l, abs(m)]
+        elseif m > 0
+            a_matrix[m, -m] = alm_coeff[l, abs(m)]
+            a_matrix[m,  m] = alm_coeff[l, abs(m)] * (-1)^abs(m)
+        else
+            a_matrix[m, -abs(m)] = alm_coeff[l, abs(m)] * 1im
+            a_matrix[m,  abs(m)] = alm_coeff[l, abs(m)] * (-1)^m * -1im
+        end
+    end
+    return parent(a_matrix)
+end

test/runtests.jl

@@ -1,6 +1,106 @@
 using CEF
-using Test
+using LinearAlgebra, Random, Test
 
-@testset "CEF.jl" begin
-    # Write your tests here.
+const PREC::Float64 = 1.0e-7
+const SDIG::Int64 = 7
+
+"""
+    test_mcphase()
+
+Returns true if extended Stevens operators for up to J=7/2 are equal to the
+operators tabulated in the mcphase manual up to a numerical precision set by PREC.
+"""
+function test_mcphase()::Bool
+    ions = ["Ce3+", "Pr3+", "Nd3+", "Pm3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+"]
+    for ion in ions
+        for l in [2, 4, 6]
+            for m in -l:1:l
+                @assert isapprox(CEF.stevens_EO(ion, l, m), CEF.stevens_O(ion, l, m), atol=PREC)
+            end
+        end
+    end
+    return true
+end
+
+
+"""
+    wignerD_unitary()
+
+Test whether the Wigner D matrix is unitary for Euler angles alpha and beta
+in 0, 2pi an 0, pi respectively.
+Values of l in [2, 4, 6] are tested since are most relevant for spectroscopy.
+"""
+function wignerD_unitary()::Bool
+    for l in [2, 4, 6]
+        for alpha in 0:0.05:2pi
+            for beta in 0:0.05:pi
+                @assert isunitary(CEFOracle.wigner_D(l, alpha, beta, 0))
+            end
+        end
+    end
+    return true
+end
+
+
+"""
+    rotation_invariance_eigenvalues()
+
+Upon rotation of the CEF Hamiltonian, the calculated eigenvalues (and
+eigenvectors up to a phase factor) should remain invariant.
+This assertion is made for a particular case of CEF Hamiltonian rotated
+by Euler angles alpha and beta in 0, 2pi and 0, pi respectively
+"""
+function rotation_invariance_eigenvalues(; seed=777)::Bool
+    rng = MersenneTwister(seed)
+    ions = ["Ce3+", "Pr3+", "Nd3+", "Pm3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+"]
+    bfactors = blm_dframe(Dict(
+                    "B2m2"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B2m1"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B20"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B21"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B22"=>     rand(rng, -0.1:1e-7:0.1),
+
+                    "B4m4"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B4m3"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B4m2"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B4m1"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B40"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B41"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B42"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B43"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B44"=>     rand(rng, -0.1:1e-7:0.1),
+
+                    "B6m6"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B6m5"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B6m4"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B6m3"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B6m2"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B6m1"=>    rand(rng, -0.1:1e-7:0.1),
+                    "B60"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B61"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B62"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B63"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B64"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B65"=>     rand(rng, -0.1:1e-7:0.1),
+                    "B66"=>     rand(rng, -0.1:1e-7:0.1),
+    ))
+    for ion in ions
+        mag_ion = single_ion(ion)
+        evals_0 = eigvals(cef_hamiltonian(mag_ion, bfactors)) # original eigenvalues
+        for alpha in 0:0.05:2pi
+            for beta in 0:0.05:pi
+                # calculate eigenvalues of system in rotated coordinate frame
+                # they should be the same as in the non-rotated system
+                @assert isapprox(eigvals(cef_hamiltonian(mag_ion, rotate_blm(bfactors, alpha, beta, 0))), evals_0, atol=PREC)
+            end
+        end
+    end
+    return true
 end
+
+
+@testset "CEF.jl" begin
+    @test test_mcphase()
+    @test wignerD_unitary()
+    @test rotation_invariance_eigenvalues()
+end