Different settings of the Oh point group appear incorrect

Dear Ted Trewick,

thanks for your comment. I'm not able to reproduce your error. Can you sent me your Quanty script where things go wrong.

Indeed as you noticed the different settings refer to rotations of the crystal field or point group. Within a crystal there are often many different orientations possible, especially if one has a cubic crystal with small local distortions. This also means the eigenvalues should be independent of the orientation. In the help files you can see, for each setting the directions of the high symmetry axis. In your calculations you need to look for the different Wickoff positions and determine the local symmetry axis for each position and then relate them to the orientation you have.

Note that the 111z direction has the C4 in the 111 direction and the C3 in the z direction. This orientation has no high symmetry axis in the y direction and should be complex.

Below you find a script that iterates over a few settings and calculates the eigenvalues. This looks o.k. for the documented and implemented functions in Quanty.

Best wishes, Maurits

And here the script to test

-- This scripts compares the different settings for point-groups as defined on the Quanty web-page.

print("Oh")
print("We take \"Random\" values for the crystal fields. Note that these have no physical meaning")
Et1u = 1.0
Et2u = 0.345
Ea2u = 0.186

Akm={}
Akm.XYZ={{0, 0, (1/7)*(Ea2u + (3)*(Et1u + Et2u))} ,
	{4, 0, (-3/4)*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)} ,
 {4,-4, (-3/4)*((sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {4, 4, (-3/4)*((sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {6, 0, (39/280)*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u))} ,
 {6,-4, (-39/40)*((1/(sqrt(14)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 4, (-39/40)*((1/(sqrt(14)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} }

Akm.OneOneOnez={{0, 0, (1/7)*(Ea2u + (3)*(Et1u + Et2u))} ,
	{4, 0, (1/2)*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)} ,
 {4, 3, (-1/2+-1/2*I)*((sqrt(5/7))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {4,-3, (1/2+-1/2*I)*((sqrt(5/7))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {6, 0, (26/105)*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u))} ,
 {6,-3, (-13/12+13/12*I)*((1/(sqrt(105)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 3, (13/12+13/12*I)*((1/(sqrt(105)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6,-6, (-13/60*I)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 6, (13/60*I)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} }

Akm.ZeroSqrt2minu1z = {{0, 0, (1/7)*(Ea2u + (3)*(Et1u + Et2u))} ,
	{4, 0, (1/2)*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)} ,
 {4,-3, (-I)*((sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {4, 3, (-I)*((sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {6, 0, (26/105)*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u))} ,
 {6,-3, (13/6*I)*((1/(sqrt(210)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 3, (13/6*I)*((1/(sqrt(210)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6,-6, (-13/60)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 6, (-13/60)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} }

Akm.sqrt201z = {{0, 0, (1/7)*(Ea2u + (3)*(Et1u + Et2u))} ,
	{4, 0, (1/2)*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)} ,
 {4,-3, (-1)*((sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u))} ,
 {4, 3, (sqrt(5/14))*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)} ,
 {6, 0, (26/105)*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u))} ,
 {6, 3, (-13/6)*((1/(sqrt(210)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6,-3, (13/6)*((1/(sqrt(210)))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6,-6, (13/60)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} ,
 {6, 6, (13/60)*((sqrt(11/21))*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u)))} }


print(Akm)

NF = 14
NB = 0
IndexDn_4f = {0,2,4,6,8,10,12}
IndexUp_4f = {1,3,5,7,9,11,13}
H={}
for k, v in pairs(Akm) do
	H[k] = NewOperator("CF",NF,IndexDn_4f,IndexUp_4f,v)
end

print(H)

val={}
fun={}
for k, v in pairs(H) do
	val[k], fun[k] = Eigensystem(OperatorToMatrix(v))
end

print(val)

Dear Maurits, Many thanks for your reply.

Here is a MWE for my issue, modified from yours, also with potentials taken from the Quanty docs. I have been using the other form of the Akm parameters as in my use case they are to be estimated from a point charge model.

Thanks, Ted

Here is the code to calculate the eigenvalues, which I expect to be equal

local A = {
    [0] = 0,
    [4] = 0.2,
    [6] = 0.008
}

AkmOh={}

AkmOh.XYZ = {{0, 0, A[0]} , 
       {4, 0, A[4]} , 
       {4,-4, (sqrt(5/14))*(A[4])} , 
       {4, 4, (sqrt(5/14))*(A[4])} , 
       {6, 0, A[6]} , 
       {6,-4, (-1)*((sqrt(7/2))*(A[6]))} , 
       {6, 4, (-1)*((sqrt(7/2))*(A[6]))} }
AkmOh.oneoneonez = {{0, 0, A[0]} , 
       {4, 0, A[4]} , 
       {4, 3, (-1+-1*I)*((sqrt(5/7))*(A[4]))} , 
       {4,-3, (1+-1*I)*((sqrt(5/7))*(A[4]))} , 
       {6, 0, A[6]} , 
       {6,-3, (-1/8+1/8*I)*((sqrt(35/3))*(A[6]))} , 
       {6, 3, (1/8+1/8*I)*((sqrt(35/3))*(A[6]))} , 
       {6,-6, (-1/8*I)*((sqrt(77/3))*(A[6]))} , 
       {6, 6, (1/8*I)*((sqrt(77/3))*(A[6]))} }

NF=14
NB=0
IndexDn_4f={0,2,4,6,8,10,12}
IndexUp_4f={1,3,5,7,9,11,13}

H={}
for k, v in pairs(AkmOh) do
	H[k] = NewOperator("CF",NF,IndexDn_4f,IndexUp_4f,v)
end


val={}
fun={}
for k, v in pairs(H) do
	val[k], fun[k] = Eigensystem(OperatorToMatrix(v))
end

print(val)

and here some python code to plot the potential in question

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from scipy.special import sph_harm_y
from numpy import sqrt

A0 = 1.0
A4 = 0.8
A6 = 0.05

Akm = {}

Akm['XYZ'] = [
    [0, 0, A0] , 
    [4, 0, A4] , 
    [4,-4, (np.sqrt(5/14))*(A4)] , 
    [4, 4, (np.sqrt(5/14))*(A4)] , 
    [6, 0, A6] , 
    [6,-4, (-1)*((np.sqrt(7/2))*(A6))] , 
    [6, 4, (-1)*((np.sqrt(7/2))*(A6))] ]


Akm['111z_Quanty'] = [ # 111z
    (0,  0, A0),
    (4,  0, A4),
    (4,  3, (-1 - 1j) * (np.sqrt(5/7) * A4)),
    (4, -3, ( 1 - 1j) * (np.sqrt(5/7) * A4)),
    (6,  0, A6),
    (6, -3, (-1/8 + 1j/8) * (np.sqrt(35/3) * A6)),
    (6,  3, ( 1/8 + 1j/8) * (np.sqrt(35/3) * A6)),
    (6, -6, (-1j/8) * (np.sqrt(77/3) * A6)),
    (6,  6, (1j/8) * (np.sqrt(77/3) * A6)) ]


Akm['111z_Ted'] = [
        [0,  0, A0],
        [4, -3,  2/21*sqrt(70)*(-1)**(3/4)*A4],
        [4,  0,                      -2/3*A4],
        [4,  3,  2/21*sqrt(70)*(-1)**(1/4)*A4],
        [6, -6,         -2/27*1j*sqrt(231)*A6],
        [6, -3, 2/27*sqrt(210)*(-1)**(3/4)*A6],
        [6,  0,                      16/9*A6],
        [6,  3, 2/27*sqrt(210)*(-1)**(1/4)*A6],
        [6,  6,          2/27*1j*sqrt(231)*A6],
        ]



fig, axes = plt.subplots(figsize=(14,5), ncols=3, subplot_kw={'projection': '3d'})


theta = np.linspace(0, np.pi, 50)     # polar angle
phi = np.linspace(0, 2 * np.pi, 50)   # azimuthal angle
theta, phi = np.meshgrid(theta, phi)


for ax, akm_name in zip(axes, Akm):
    y_total = np.zeros_like(theta, dtype=complex)
    for l, m, coeff in Akm[akm_name]:
        y_total += coeff * np.sqrt(4*np.pi/(2*l+1)) * sph_harm_y(l, m, theta, phi)
    
    r = y_total.real
    x = r * np.sin(theta) * np.cos(phi)
    y = r * np.sin(theta) * np.sin(phi)
    z = r * np.cos(theta)
    
    surf = ax.plot_surface(x, y, z,
                           rstride=1, cstride=1, linewidth=0, antialiased=False)
    
    ax.set_title(akm_name)
    ax.set_aspect('equal')
plt.show()

Dear Ted,

I might see where your problem comes from. For each orientation and pointgroup I write the potential as $$V(r,\theta,\phi) = \sum_{k=0}^{\infty} \sum_{m=-k}^{k} A_{k,m}(r) C^{(m)}_k(\theta,\phi)$$ which also means that for each orientation the physical meaning of these expansion coefficients is different. In Oh point group there are relations between the expansion coefficients, again different for each setting. To make things short, the A[0], A[4], and A[6] in the XYZ setting are not the same parameters as the A[0], A[4] and A[6] in the 111z setting. The second part of the help files relates the expansion expanded on spherical harmonics to the eigen energies of the point group.

For the plotting, you probably run into the issue that for degenerate orbitals (the t1u and t2u representations) any linear combination of these orbitals is a fine solution. Such that the orbitals can look very different although they span the same space.

To be complete, here my changes to your script to relate the Akm's for an Oh expansion with the high symmetry axis in the X, Y and Z direction to an expansion with the C4 axis in the 111 direction and the C3 in the z direction

PS you can find all relations on the pages https://www.quanty.org/physics_chemistry/point_groups/oh/orientation_111z and https://www.quanty.org/physics_chemistry/point_groups/oh/orientation_xyz

especially if you expand the sections below Potential for f orbitals, or if you look at the one particle Hamiltonian matrix on a basis of cubic harmonics

local A = {
		[0] = 0,
		[4] = 0.2,
		[6] = 0.008
}

AkmOh={}

AkmOh.XYZ = {{0, 0, A[0]} ,
			 {4, 0, A[4]} ,
			 {4,-4, (sqrt(5/14))*(A[4])} ,
			 {4, 4, (sqrt(5/14))*(A[4])} ,
			 {6, 0, A[6]} ,
			 {6,-4, (-1)*((sqrt(7/2))*(A[6]))} ,
			 {6, 4, (-1)*((sqrt(7/2))*(A[6]))} }


print("Eigen energies from paramters in XYZ orientation ")

Ea2u = A[0] - (12/33) * A[4] + (240/429)*A[6]
Et1u = A[0] + ( 6/33) * A[4] + (100/429)*A[6]
Et2u = A[0] - ( 2/33) * A[4] - (180/429)*A[6]

print("Ea2u ",Ea2u)
print("Et1u ",Et1u)
print("Et2u ",Et2u)

local A = {
		[0] = 0,
		[4] = 0.2,
		[6] = 0.008
}
print("Potential expanded on spherical harmonics for the 111 z direction C4 // 111 -- C3 // z")

A[0] = (1/7)*(Ea2u + (3)*(Et1u + Et2u))
A[4] = (1/2)*((2)*(Ea2u) + (-3)*(Et1u) + Et2u)
A[6] = (26/105)*((4)*(Ea2u) + (5)*(Et1u) + (-9)*(Et2u))

print("A[0] ", A[0])
print("A[4] ", A[4])
print("A[6] ", A[6])
			
AkmOh.oneoneonez = {{0, 0, A[0]} ,
			 {4, 0, A[4]} ,
			 {4, 3, (-1+-1*I)*((sqrt(5/7))*(A[4]))} ,
			 {4,-3, (1+-1*I)*((sqrt(5/7))*(A[4]))} ,
			 {6, 0, A[6]} ,
			 {6,-3, (-1/8+1/8*I)*((sqrt(35/3))*(A[6]))} ,
			 {6, 3, (1/8+1/8*I)*((sqrt(35/3))*(A[6]))} ,
			 {6,-6, (-1/8*I)*((sqrt(77/3))*(A[6]))} ,
			 {6, 6, (1/8*I)*((sqrt(77/3))*(A[6]))} }

NF=14
NB=0
IndexDn_4f={0,2,4,6,8,10,12}
IndexUp_4f={1,3,5,7,9,11,13}

H={}
for k, v in pairs(AkmOh) do
	H[k] = NewOperator("CF",NF,IndexDn_4f,IndexUp_4f,v)
end


val={}
fun={}
for k, v in pairs(H) do
	val[k], fun[k] = Eigensystem(OperatorToMatrix(v))
end

print(val)

Dear Maurits,

You are of course correct; I now see this was totally the result of a misconception on my part. Thank you so much for clarifying the meaning of the A[i] parameters.

I think for now I will continue to find the rotated potentials the old fashioned way, using the Wigner-D matrices, mostly to avoid the mistakes I'm sure to make when converting by expanding the potentials on the representations.

Many thanks,

Ted

Dear Ted,

No problem. Be careful if you rotate potentials and or Hamiltonians though. Do you rotate the basis or the vector, i.e. do you need $U^{\dagger}$, $U^{T}$, or $U^{*}$ as your rotation matrix.

PS you can test in Quanty with the function Rotate, which can rotate Wavefunctions, operators and matrices. https://www.quanty.org/documentation/language_reference/functions/rotate

PS2 the reason I have the documentation is because I made so many errors on this that I started to document and test the different rotations and orientations of the different crystal fields. The pages that are implemented relate to point groups of materials I once studied….

Best wishes, Maurits