Verify that anisotropy operators rotate correctly

[kbarros]

David and I encountered this surprising Sunny behavior. I need to verify that it is doing the correct thing.

R = Sunny.rotation_between_vectors([1,1,1], [0,0,1])
𝒪 = stevens_matrices(Inf)

latvecs = hcat([0,1/2,1/2], [1/2, 0, 1/2], [1/2,1/2,0])
cryst1 = Crystal(latvecs, [[0,0,0]], 166)
view_crystal(cryst1)
print_site(cryst1, 1)

# At quartic level:
op1 = (-7𝒪[4,-3]-2𝒪[4,-2]+𝒪[4,-1]+𝒪[4,1]+7𝒪[4,3]) / (14/3)
op2 = (𝒪[4,0]+5𝒪[4,4]) / (-40/3)

cryst2 = Crystal(R*latvecs, [[0,0,0]], 166)
view_crystal(cryst2)
print_site(cryst2, 1)
print_site(cryst1, 1; R) # equivalent

# At quartic level:
op3 = 𝒪[4,0] * 9
op4 = (𝒪[4,-3]-𝒪[4,3]) * 18

# Here I expected to recreate op3 and op4, but the result is different
print_stevens_expansion(rotate_operator(op1, R))
print_stevens_expansion(rotate_operator(op2, R))

[kbarros]

It is correct. With the appropriate linear combination of rotated op1 and op2, one can reproduce op3 and op4 above.

using Sunny
R = Sunny.rotation_between_vectors([1,1,1], [0,0,1])
𝒪 = stevens_matrices(Inf)

latvecs = hcat([0,1/2,1/2], [1/2, 0, 1/2], [1/2,1/2,0])
cryst1 = Crystal(latvecs, [[0,0,0]], 166)
print_site(cryst1, 1)
op1 = (-7𝒪[4,-3]-2𝒪[4,-2]+𝒪[4,-1]+𝒪[4,1]+7𝒪[4,3]) / (14/3)
op2 = (𝒪[4,0]+5𝒪[4,4]) / (-40/3)

op1_rot = rotate_operator(op1, R')
op2_rot = rotate_operator(op2, R')
print_stevens_expansion((op2_rot - op1_rot)*(140/27)) # 𝒪₄₀
print_stevens_expansion(op1_rot + op2_rot - (13/27) * (op1_rot - op2_rot)) # 2𝒪₄₋₃ - 2𝒪₄₃

# Can get the same operators by performing symmetry analysis on the rotated
# crystal
cryst2 = Crystal(R*latvecs, [[0,0,0]], 166)
print_site(cryst2, 1)
print_site(cryst1, 1; R) # equivalent