Directional X-ray Spectroscopy, Potential Bug
asked by Charles Cardot (2023/04/07 23:56)
Hello,
I have been experimenting with calculating direction dependent x-ray emission using Quanty and the DFT+MLFT framework. However, I seem to be running into an issue with the mean field operators calculated from the Tight Binding (TB) object.
The directional dependence is encoded in the DFT step, where I distort the crystal lattice slightly (~2%) along a particular axis.
# lattice constants; @l@ 4.27633759 4.18633759 4.18633759 # set axis angles @a@ 90. 90. 90. # setup Wyckoff positions @n@2 # Now, give list of ALL !!! Wyckoff positions. @1@ Ni @ 0.00000000 0.00000000 0.00000000 @2@ O @ 0.00000000 0.00000000 0.50000000
When I stretch along the a-axis and go through the regular Wannierization process to create a tight binding Hamiltonian, I get an HDFT that looks like
a_stretched_HDFT: https://drive.google.com/file/d/1sE72v2rbbPM9Bof4Hg_CgNgvqmtjKjYY/view?usp=share_link
which makes sense in the context of the distorted octahedral crystal field. The xy and xz orbitals are still degenerate, while the yz is now lower energy. Also the z^2-r^2 and x^2-y^2 orbitals become mixed. The MLFT calculation within Quanty proceeds much the same way as described in the 2022 and 2019 Heidelberg tutorials for DFT + MLFT, where the mean field components Coulomb interaction is subtracted from HDFT to avoid double counting. The full code that I am running can be found here (https://drive.google.com/drive/folders/1G0iSTfIMT8i39CWWg1StEs5Kj26ZYQ1W?usp=share_link), but an excerpt is shown below.
print("--Define Intermediate State Hamiltonian--\n")
--Hamiltonian = HDFT + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
Hamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
Hamiltonian = Hamiltonian + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
Hamiltonian = Hamiltonian + zeta_3d * Oppldots_3d
Hamiltonian = Hamiltonian + F0sd*OppF0sd
Hamiltonian = Hamiltonian/2
Hamiltonian = Hamiltonian + ConjugateTranspose(Hamiltonian)
esinter = -nd*Usd
edinter = -((-20*Delta + 19*nd*Udd+nd*nd*Udd+40*Usd)/(2*(10+nd)))
eLinter = nd*((1+nd)*Udd/2 - Delta + 2*Usd)/(10+nd)
OperatorSetTrace(Hamiltonian,esinter,Index["Ni_1s"])
OperatorSetTrace(Hamiltonian,edinter,Index["Ni_3d"])
OperatorSetTrace(Hamiltonian,eLinter,Index["Ligand_d"])
print("--Define XES-Hamiltonian--\n")
--XESHamiltonian = HDFT + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
XESHamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
XESHamiltonian = XESHamiltonian + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
XESHamiltonian = XESHamiltonian + zeta_3d * Oppldots_3d
XESHamiltonian = XESHamiltonian + zeta_2p * Oppldots_2p
XESHamiltonian = XESHamiltonian + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3
XESHamiltonian = XESHamiltonian/2
XESHamiltonian = XESHamiltonian + ConjugateTranspose(XESHamiltonian)
epfinal = -nd*Upd
edfinal = -((-20*Delta + 19*nd*Udd + nd*nd*Udd + 120*Upd)/(2*(10 + nd)))
eLfinal = nd*(-2*Delta + Udd + nd*Udd + 12*Upd)/(2*(10 + nd))
OperatorSetTrace(XESHamiltonian,epfinal,Index["Ni_2p"])
OperatorSetTrace(XESHamiltonian,edfinal,Index["Ni_3d"])
OperatorSetTrace(XESHamiltonian,eLfinal,Index["Ligand_d"])
print("--Compute eigenstates--\n")
-- we now can create the lowest Npsi eigenstates:
Npsi=6
StartRestrictions = {NFermi, 0, {DeterminantString(NFermi,Index["Ni_1s"]),1,1},
{DeterminantString(NFermi,Index["Ni_3d"],Index["Ligand"]),16+nd,16+nd},
{DeterminantString(NFermi,Index["Ni_2p"]),6,6}}
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi, {{'Zero',1e-12},{'Epsilon',1e-12}})
psiList = Chop(psiList)
print(StartRestrictions)
Npsiline = {}
for i=1,Npsi,1
do
table.insert(Npsiline,1)
end
print("Npsiline")
print(Npsiline)
print("--Create the Spectra--\n")
Emin = -50
Emax = 50
NE= 50000
-- Constant Lorentzian Broadening --
------------------------------------
Gamma = 0.1
------------------------------------
----------- For Sticks -------------
Spectra_z = CreateSpectra(XESHamiltonian, TXESzdag, psiList, {{"Emin",Emin}, {"Emax",Emax}, {"NE",NE}, {"Gamma",Gamma}})
Spectra_x = CreateSpectra(XESHamiltonian, TXESxdag, psiList, {{"Emin",Emin}, {"Emax",Emax}, {"NE",NE}, {"Gamma",Gamma}})
Spectra_y = CreateSpectra(XESHamiltonian, TXESydag, psiList, {{"Emin",Emin}, {"Emax",Emax}, {"NE",NE}, {"Gamma",Gamma}})
Spectra_z.Broaden(0, 1.9)
Spectra_x.Broaden(0, 1.9)
Spectra_y.Broaden(0, 1.9)
Spectra_z = Spectra.Sum(Spectra_z,Npsiline)
Spectra_x = Spectra.Sum(Spectra_x,Npsiline)
Spectra_y = Spectra.Sum(Spectra_y,Npsiline)
Spectra_z.Print({{"file", "XES_zpol.dat"}})
Spectra_x.Print({{"file", "XES_xpol.dat"}})
Spectra_y.Print({{"file", "XES_ypol.dat"}})
print("Finished")
The plotted output looks like
XES_including_MFDFT: https://drive.google.com/file/d/1hYutZdeujtK8lohFPPD7sZ-7z5EN4asl/view?usp=share_link
where clearly the x, y, and z polarized emissions are different. Given that the distortion is along the a-axis, I expect the y and z spectra to be identical, with the x spectra being slightly different. In investigating this, I noticed that omitting the subtraction of the Mean Field DFT contribution,
print("--Define Intermediate State Hamiltonian--\n")
Hamiltonian = HDFT + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
--Hamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
--Hamiltonian = Hamiltonian + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
Hamiltonian = Hamiltonian + zeta_3d * Oppldots_3d
Hamiltonian = Hamiltonian + F0sd*OppF0sd
Hamiltonian = Hamiltonian/2
Hamiltonian = Hamiltonian + ConjugateTranspose(Hamiltonian)
esinter = -nd*Usd
edinter = -((-20*Delta + 19*nd*Udd+nd*nd*Udd+40*Usd)/(2*(10+nd)))
eLinter = nd*((1+nd)*Udd/2 - Delta + 2*Usd)/(10+nd)
OperatorSetTrace(Hamiltonian,esinter,Index["Ni_1s"])
OperatorSetTrace(Hamiltonian,edinter,Index["Ni_3d"])
OperatorSetTrace(Hamiltonian,eLinter,Index["Ligand_d"])
print("--Define XES-Hamiltonian--\n")
XESHamiltonian = HDFT + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
--XESHamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
--XESHamiltonian = XESHamiltonian + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
XESHamiltonian = XESHamiltonian + zeta_3d * Oppldots_3d
XESHamiltonian = XESHamiltonian + zeta_2p * Oppldots_2p
XESHamiltonian = XESHamiltonian + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3
XESHamiltonian = XESHamiltonian/2
XESHamiltonian = XESHamiltonian + ConjugateTranspose(XESHamiltonian)
epfinal = -nd*Upd
edfinal = -((-20*Delta + 19*nd*Udd + nd*nd*Udd + 120*Upd)/(2*(10 + nd)))
eLfinal = nd*(-2*Delta + Udd + nd*Udd + 12*Upd)/(2*(10 + nd))
OperatorSetTrace(XESHamiltonian,epfinal,Index["Ni_2p"])
OperatorSetTrace(XESHamiltonian,edfinal,Index["Ni_3d"])
OperatorSetTrace(XESHamiltonian,eLfinal,Index["Ligand_d"])
causes the issue to go away and the spectra exhibits the expected behavior,
XES_excluding_MFDFT: https://drive.google.com/file/d/157Vyqxn2Ld20tajT0w3anxgFYQB9p8xu/view?usp=share_link,
with the y and z polarized emission lying exactly on top of each other, and the x emission being distinct from the other two directions.
What this tells me is that the crystal structure that is reflected in the tight binding Hamiltonain, HDFT, is indeed correct, and produces the expected directional dependence. However, something about the Coulomb correction is not behaving as expected. I'm not sure if this is an issue of numerical accuracy (the distortion is small) or a potential bug in the implementation of the MeanFieldOperator. Any help that can be provided would be much appreciated.
Answers
Update: After working extensively with the Quanty developers we have discovered that this is indeed an expected behavior based on how the DFT correction to the MeanFieldOperator is defined (https://www.quanty.org/documentation/language_reference/functions/meanfieldoperator). When the option “AddDFTSelfInteraction” is turned on, a term is added to the MeanFieldOperator which contributes only along the diagonal (m=m). If there are off diagonal coupling terms (for example the mixing of different tesseral harmonics such as in my distorted crystal case), then the mean field operator will not follow the same symmetry operations as the density matrix which was used to construct it. This leads to symmetry issues when the MeanFieldOperator is combined with HDFT later in the calculation.
I know that this won't be an issue for isotropic coordination environments (perfect Oh or perfect Td) and I have a hunch that this isn't an issue for isotropic spectra (averaging over x, y, and z polarizations) of any system. That being said, it will become noticeable for lower symmetry systems or for any use cases looking at polarized spectra such as XMCD.
An updated version of the MeanFieldOperator which sums over all m,n states instead of just m will give a rotationally invariant DFT correction. I hope that something like this will be included in the next public release of Quanty.