From DFT to MLFT Transition Operator for 2p X-ray Photoemission Spectra
asked by Hamza (2024/11/23 15:44)
Dear Quanty Developers,
I am currently following the tutorial from dfttolmft. While the tutorial provides an example for XAS calculations, I noticed there is no equivalent example for photoemission.
I attempted to use the transition operator defined in the standard ligand field theory, but the results were significantly different from what I expected. Below is a snippet of my code for reference:
--============================================================================
function YlmClmYlm (l1,m1,l2,m2,l3,m3)
return (-1)^m1 * sqrt((2*l1+1)*(2*l3+1)) * ThreeJSymbol({l1,0},{l2,0},{l3,0}) * ThreeJSymbol({l1,-m1},{l2,m2},{l3,m3})
end
function sqr(x)
return x*x
end
function TransitionOperator (l3, IndexDn, IndexUp, NF, epsilon_x, epsilon_y, epsilon_z, detector_theta, detector_phi)
local T = 0
local epsilon={}
epsilon[-1] = sqrt(1/2)*( epsilon_x - I*epsilon_y) / sqrt(sqr(epsilon_x) + sqr(epsilon_y) + sqr(epsilon_z))
epsilon[ 0] = epsilon_z / sqrt(sqr(epsilon_x) + sqr(epsilon_y) + sqr(epsilon_z))
epsilon[ 1] = sqrt(1/2)*(-epsilon_x - I*epsilon_y) / sqrt(sqr(epsilon_x) + sqr(epsilon_y) + sqr(epsilon_z))
-- T = sum_{m,m',m''} a_m <m|C_epsilon_m'|s> Y(lm'',theta,phi)
for m1=-1,1 do -- loop over angular momentum of the p electron
for m2=-1,1 do -- loop over angular momentum of the dipole operator
for m3=-l3,l3 do -- loop over angular momentum of the free electron
T = T + NewOperator("An", NF, {IndexDn[m1+2]}, {epsilon[m2] * YlmClmYlm(1,m1,1,m2,l3,m3) * SphericalHarmonicY(l3,m3,detector_theta,detector_phi)})
+ NewOperator("An", NF, {IndexUp[m1+2]}, {epsilon[m2] * YlmClmYlm(1,m1,1,m2,l3,m3) * SphericalHarmonicY(l3,m3,detector_theta,detector_phi)})
end
end
end
return Chop(T)
end
--============================================================================
.
.
.
.
--=========================================================
epsilon = {1,0,0}
detector_theta = pi/4
detector_phi = 0
IndexDn_2p={ 0, 2, 4}
IndexUp_2p={ 1, 3, 5}
NF = 26
TPES_s = TransitionOperator (0,IndexDn_2p, IndexUp_2p, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TPES_d = TransitionOperator (2, IndexDn_2p, IndexUp_2p, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TPES = TPES_s + TPES_d
-- calculating the spectra is simple and single line once all operators and wave-functions
-- are defined.
cPESSpectra = CreateSpectra(XASHamiltonian, {TPES_s,TPES_d,TPES}, psiList, {{"Emin",-30}, {"Emax",10}, {"NE",4000}, {"Gamma",0.1}})
Could you please confirm if I am on the right track? If not, I would greatly appreciate your guidance on how to properly define the transition operator for photoemission in this case. Additionally, are there any other changes I need to make to the 3_XAS file to accommodate this calculation?
Thank you for your support!
Best regards, Hamza
Answers
Dear Hamza,
Indeed changing the transition operator could be sufficient.
For core level XPS you seem to be on the right track. Note that the full operator you measure in the experiment should be c_1 * TPES_s + c_2 * TPES_d with c_1 and c_2 complex parameters that depend on the photon energy and core level involved (these are determined by integrals over the radial matrix elements).
I'm not sure how you defined the onsite energies of the different configurations or shells, i.e. the relation between the epsilon_d, epsilon_L and U, Delta. This is the first part I would check if there is something strange going on.
What is different between the calculation and your expectations?
Best wishes, Maurits
Dear Maurits,
Thank you for your advice.
I discovered an error: I was using the wrong index for spin-up and spin-down in the transition operator. However, after fixing it, I am pleased with the result for NiO it was similar to the results from the ligand field directory.
Regarding the onsite energy, I did not change it. Do you think I need to modify it in the case of Iron? Just as a reminder, I am using the DFT-to-LMFT method.
Best regards, Hamza
Dear Hamza,
no need to change anything in the general setup of the script. You do have the parameter n which should be changed according to the filling / element you use.
Best wishes, Maurits
Dear Maurits,
Thank you for your reply. I came across the function CreateAtomicIndicesDict. Do we need to compute the number of the ligand after the rotation of the created cluster Hamiltonian? Either for XAS or XPS.
print("Create Cluster Hamiltonian after rotation\n") HDFT, ClusterTB, IndexFPLO = CreateClusterHamiltonian(RTB, {"open", RTB.Atoms}, {{"AddSpin",true}, {"ReturnTBIndicesDict",true}}) --Now we want to add the Ni 2_p shell to the system --New index object (first argument is a list of all orbitals, the second argument groups the two Ni orbitals) Index, NFermi = CreateAtomicIndicesDict({"Ni_3d","Ligand_d","Ni_2p"}, {{"Ni",{"Ni_2p","Ni_3d"}}}) HDFT.NF = NFermi print("The indices of the system are now fixed:")Best regard, Hamza