Skip to content

Commit 8d63e76

Browse files
abehersanabehersan
authored
Merge pull request #2 from abehersan/dev
Dev
2 parents 18e6fec + 81c10ef commit 8d63e76

4 files changed

Lines changed: 252 additions & 4 deletions

File tree

‎.gitignore‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,2 @@
11
./Manifest.toml
2-
./github
2+
./.github

‎src/CEF.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,4 +65,7 @@ include("./cef_neutronxsection.jl")
6565
export cef_neutronxsection_crystal!
6666
export cef_neutronxsection_powder!
6767

68+
include("./cef_rotation.jl")
69+
export rotate_blm
70+
6871
end

‎src/cef_rotation.jl‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
function rotation_unitary(J::Real, n::Vector{<:Real}, phi::Real)::Matrix{ComplexF64}
2+
Jx = spin_operators(J, "x")
3+
Jy = spin_operators(J, "y")
4+
Jz = spin_operators(J, "z")
5+
nnorm = normalize(n)
6+
Jproj = sum(nnorm .* [Jx, Jy, Jz])
7+
U = exp(-1im * phi * Jproj)
8+
@assert isunitary(U)
9+
return U
10+
end
11+
12+
13+
function ZYZ_rotmatrix(alpha::Real, beta::Real, gamma::Real)::Matrix{Float64}
14+
a, b, g = map(Float64, [alpha, beta, gamma])
15+
return Z_rot(a) * Y_rot(b) * Z_rot(g)
16+
end
17+
18+
19+
X_rot(theta::Float64)::Matrix{Float64} = [1 0 0;
20+
0 cos(theta) -sin(theta);
21+
0 sin(theta) cos(theta)]
22+
23+
24+
Y_rot(theta::Float64)::Matrix{Float64} = [cos(theta) 0 sin(theta);
25+
0 1 0;
26+
-sin(theta) 0 cos(theta)]
27+
28+
29+
Z_rot(theta::Float64)::Matrix{Float64} = [cos(theta) -sin(theta) 0;
30+
sin(theta) cos(theta) 0;
31+
0 0 1]
32+
33+
34+
function get_euler_angles(v::Vector{<:Real})::Tuple{Float64, Float64, Float64}
35+
if isequal(zeros(Real, 3), v)
36+
return (0.0, 0.0, 0.0)
37+
end
38+
v_norm = v / norm(v)
39+
alpha = atan(v_norm[2], v_norm[1])# / pi * 180.0
40+
beta = atan((v_norm[1]^2 + v_norm[2]^2), v_norm[3])# / pi * 180.0
41+
gamma = 0.0
42+
return (alpha, beta, gamma) # in radian
43+
end
44+
45+
46+
function wigner_D(l::Int, alpha::Real, beta::Real, gamma::Real)::Matrix{ComplexF64}
47+
mdim = Int(2*l+1)
48+
ml = -l:1:l
49+
rotmat = Matrix{ComplexF64}(undef, (mdim, mdim))
50+
for (idmp, mm) in enumerate(ml)
51+
for (idm, mp) in enumerate(ml)
52+
rotmat[idmp, idm] = (-1)^(mp - mm) * rotation_matrix_element(l, mm, mp, alpha, beta, gamma)
53+
end
54+
end
55+
return round.(rotmat, digits=SDIG)
56+
end
57+
58+
59+
function rotation_matrix_element(l::Int, m::Int, mp::Int, alpha::Real, beta::Real, gamma::Real)::ComplexF64
60+
return exp(-1.0im*mp*alpha) * small_d(l, mp, m, beta) * exp(-1.0im*m*gamma)
61+
end
62+
63+
64+
function small_d(l::Int, mp::Int, m::Int, beta::Real)::Float64
65+
if iszero(beta)
66+
return isequal(m, mp) * 1.0 # delta function d_m'm if beta=0
67+
end
68+
djmpm = 0.0
69+
smin = Int(maximum([0, -(m+mp)]))
70+
smax = Int(minimum([l-m, l-mp]))
71+
for s in smin:1:smax
72+
djmpm += (-1)^(l-m-s)*
73+
((cos(beta/2))^(2s+mp+m))*
74+
((sin(beta/2))^(2l-2s-m-mp))*
75+
binomial(Int(l+m), Int(l-mp-s))*
76+
binomial(Int(l-m), Int(s))
77+
end
78+
djmpm *= sqrt((factorial(Int(l+mp))*factorial(Int(l-mp)))/
79+
(factorial(Int(l+m))*factorial(Int(l-m))))
80+
if isequal(djmpm, NaN) || isequal(djmpm, Inf)
81+
err_message =
82+
"Matrix element djmpm is NaN or Inf!\n"*
83+
"Likely there's an issue in the inputted values of l, m or mp."
84+
@error err_message
85+
else
86+
return djmpm
87+
end
88+
end
89+
90+
91+
function rotate_blm(cefparams::DataFrame, alpha::Real, beta::Real, gamma::Real)::DataFrame
92+
cefparams_rotated = DataFrame(B = Float64[], l = Int[], m = Int[])
93+
ls = sort(collect(Set(cefparams[:, :l])))
94+
for l in ls
95+
S_matrix = rotate_stevens(l, alpha, beta, gamma)
96+
cefparams_ori = zeros(Float64, Int(2l+1)) # original CEF parameters
97+
cefparams_rot = zeros(ComplexF64, Int(2l+1)) # rotated CEF params (complex)
98+
cefparams_res = zeros(Float64, Int(2l+1)) # rotated CEF parameters (real)
99+
for (i, m) in enumerate(-l:1:l)
100+
try
101+
cefparams_ori[i] = cefparams[(cefparams.l .== l) .& (cefparams.m .== m), :B][1]
102+
catch ex
103+
if isa(ex, BoundsError)
104+
cefparams_ori[i] = 0.0
105+
end
106+
end
107+
end
108+
# S * cefparams where cefparams is a (2l+1) vector
109+
cefparams_rot .= S_matrix * cefparams_ori
110+
@assert norm(imag(cefparams_rot)) < PREC
111+
cefparams_res .= real(cefparams_rot)
112+
append!(cefparams_rotated, DataFrame("B"=>cefparams_res, "l"=>fill(l, Int(2l+1)), "m"=>-l:1:l))
113+
end
114+
return cefparams_rotated
115+
end
116+
117+
118+
function rotate_stevens(l::Int, alpha::Real, beta::Real, gamma::Real)::Matrix{ComplexF64}
119+
return transpose(inv(Alm_matrix(l))) * wigner_D(l, alpha, beta, gamma)' * transpose(Alm_matrix(l))
120+
end
121+
122+
123+
function Alm_matrix(l::Int)::Matrix{ComplexF64}
124+
alm_coeff = OffsetArray(SMatrix{6, 7}([
125+
1 1/sqrt(2) 0 0 0 0 0;
126+
sqrt(6) 1/2 1 0 0 0 0;
127+
sqrt(10) sqrt(10/3) 1/sqrt(3) sqrt(2) 0 0 0;
128+
2*sqrt(70) sqrt(7/2) sqrt(7) 1/sqrt(2) 2 2 0;
129+
6*sqrt(14) 2*sqrt(21/5) sqrt(3/5) 6*sqrt(2/5) 2/sqrt(5) 2*sqrt(2) 0;
130+
4*sqrt(231) sqrt(22) 4*sqrt(11/5) 2*sqrt(11/5) 4*sqrt(11/6) 2/sqrt(3) 4;
131+
]), 1:6, 0:6)
132+
a_matrix = OffsetArray(zeros(ComplexF64, (2l+1, 2l+1)), -l:l, -l:l)
133+
for m in -l:1:l
134+
if iszero(m)
135+
a_matrix[m, m] = alm_coeff[l, abs(m)]
136+
elseif m > 0
137+
a_matrix[m, -m] = alm_coeff[l, abs(m)]
138+
a_matrix[m, m] = alm_coeff[l, abs(m)] * (-1)^abs(m)
139+
else
140+
a_matrix[m, -abs(m)] = alm_coeff[l, abs(m)] * 1im
141+
a_matrix[m, abs(m)] = alm_coeff[l, abs(m)] * (-1)^m * -1im
142+
end
143+
end
144+
return parent(a_matrix)
145+
end

‎test/runtests.jl‎

Lines changed: 103 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,106 @@
11
using CEF
2-
using Test
2+
using LinearAlgebra, Random, Test
33

4-
@testset "CEF.jl" begin
5-
# Write your tests here.
4+
const PREC::Float64 = 1.0e-7
5+
const SDIG::Int64 = 7
6+
7+
"""
8+
test_mcphase()
9+
10+
Returns true if extended Stevens operators for up to J=7/2 are equal to the
11+
operators tabulated in the mcphase manual up to a numerical precision set by PREC.
12+
"""
13+
function test_mcphase()::Bool
14+
ions = ["Ce3+", "Pr3+", "Nd3+", "Pm3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+"]
15+
for ion in ions
16+
for l in [2, 4, 6]
17+
for m in -l:1:l
18+
@assert isapprox(CEF.stevens_EO(ion, l, m), CEF.stevens_O(ion, l, m), atol=PREC)
19+
end
20+
end
21+
end
22+
return true
23+
end
24+
25+
26+
"""
27+
wignerD_unitary()
28+
29+
Test whether the Wigner D matrix is unitary for Euler angles alpha and beta
30+
in 0, 2pi an 0, pi respectively.
31+
Values of l in [2, 4, 6] are tested since are most relevant for spectroscopy.
32+
"""
33+
function wignerD_unitary()::Bool
34+
for l in [2, 4, 6]
35+
for alpha in 0:0.05:2pi
36+
for beta in 0:0.05:pi
37+
@assert isunitary(CEFOracle.wigner_D(l, alpha, beta, 0))
38+
end
39+
end
40+
end
41+
return true
42+
end
43+
44+
45+
"""
46+
rotation_invariance_eigenvalues()
47+
48+
Upon rotation of the CEF Hamiltonian, the calculated eigenvalues (and
49+
eigenvectors up to a phase factor) should remain invariant.
50+
This assertion is made for a particular case of CEF Hamiltonian rotated
51+
by Euler angles alpha and beta in 0, 2pi and 0, pi respectively
52+
"""
53+
function rotation_invariance_eigenvalues(; seed=777)::Bool
54+
rng = MersenneTwister(seed)
55+
ions = ["Ce3+", "Pr3+", "Nd3+", "Pm3+", "Sm3+", "Tb3+", "Dy3+", "Ho3+", "Er3+", "Tm3+", "Yb3+"]
56+
bfactors = blm_dframe(Dict(
57+
"B2m2"=> rand(rng, -0.1:1e-7:0.1),
58+
"B2m1"=> rand(rng, -0.1:1e-7:0.1),
59+
"B20"=> rand(rng, -0.1:1e-7:0.1),
60+
"B21"=> rand(rng, -0.1:1e-7:0.1),
61+
"B22"=> rand(rng, -0.1:1e-7:0.1),
62+
63+
"B4m4"=> rand(rng, -0.1:1e-7:0.1),
64+
"B4m3"=> rand(rng, -0.1:1e-7:0.1),
65+
"B4m2"=> rand(rng, -0.1:1e-7:0.1),
66+
"B4m1"=> rand(rng, -0.1:1e-7:0.1),
67+
"B40"=> rand(rng, -0.1:1e-7:0.1),
68+
"B41"=> rand(rng, -0.1:1e-7:0.1),
69+
"B42"=> rand(rng, -0.1:1e-7:0.1),
70+
"B43"=> rand(rng, -0.1:1e-7:0.1),
71+
"B44"=> rand(rng, -0.1:1e-7:0.1),
72+
73+
"B6m6"=> rand(rng, -0.1:1e-7:0.1),
74+
"B6m5"=> rand(rng, -0.1:1e-7:0.1),
75+
"B6m4"=> rand(rng, -0.1:1e-7:0.1),
76+
"B6m3"=> rand(rng, -0.1:1e-7:0.1),
77+
"B6m2"=> rand(rng, -0.1:1e-7:0.1),
78+
"B6m1"=> rand(rng, -0.1:1e-7:0.1),
79+
"B60"=> rand(rng, -0.1:1e-7:0.1),
80+
"B61"=> rand(rng, -0.1:1e-7:0.1),
81+
"B62"=> rand(rng, -0.1:1e-7:0.1),
82+
"B63"=> rand(rng, -0.1:1e-7:0.1),
83+
"B64"=> rand(rng, -0.1:1e-7:0.1),
84+
"B65"=> rand(rng, -0.1:1e-7:0.1),
85+
"B66"=> rand(rng, -0.1:1e-7:0.1),
86+
))
87+
for ion in ions
88+
mag_ion = single_ion(ion)
89+
evals_0 = eigvals(cef_hamiltonian(mag_ion, bfactors)) # original eigenvalues
90+
for alpha in 0:0.05:2pi
91+
for beta in 0:0.05:pi
92+
# calculate eigenvalues of system in rotated coordinate frame
93+
# they should be the same as in the non-rotated system
94+
@assert isapprox(eigvals(cef_hamiltonian(mag_ion, rotate_blm(bfactors, alpha, beta, 0))), evals_0, atol=PREC)
95+
end
96+
end
97+
end
98+
return true
699
end
100+
101+
102+
@testset "CEF.jl" begin
103+
@test test_mcphase()
104+
@test wignerD_unitary()
105+
@test rotation_invariance_eigenvalues()
106+
end

0 commit comments

Comments
 (0)