-- This tutorial calculates the 2p to 1s x-ray emission spectra of Ni in NiO using
-- crystal field theory. We assume that the decay of the 1s core hole is fast such
-- that the 3d state is unchanged by the 1s core hole.

-- For an experiment see for example here
-- https://www.researchgate.net/figure/Ni-Ka-spectra-of-NiF2-and-K2NiF6-top-and-NiO-and-Ni2O3-bottom-measured-with_fig1_7705996

-- Within crystal-field theory the solid is approximated by a single atom in an effective 
-- potential. Although an extremely crude approximation it is useful for some cases.
-- For correlated transition metal insulators it captures the right symmetry of the 
-- localized open d-shell. It is useful to determine magnetic g-factors, energies of d-d
-- excitations or core level x-ray absorption. (2p to 3d excitations L23 edges) 

-- One should notice that the effective crystal-field potential is an affective potential
-- it is there to mimic the interaction with neighboring ligand atoms. In real materials
-- there do not exist such large electro static potentials.

-- In order to do crystal-field theory for NiO we need to define a Ni d-shell.
-- A d-shell has 10 elements and we label again the even spin-orbitals to be spin down
-- and the odd spin-orbitals to be spin up. Besides the Ni 3d shell we also need the Ni 1s
-- shell and the Ni 2p shell
NF=18
NB=0
IndexDn_1s={0}
IndexUp_1s={1}
IndexDn_2p={2,4,6}
IndexUp_2p={3,5,7}
IndexDn_3d={8,10,12,14,16}
IndexUp_3d={9,11,13,15,17}

-- just like in the previous example we define several operators acting on the Ni -3d shell

OppSx   =NewOperator("Sx"   ,NF, IndexUp_3d, IndexDn_3d)
OppSy   =NewOperator("Sy"   ,NF, IndexUp_3d, IndexDn_3d)
OppSz   =NewOperator("Sz"   ,NF, IndexUp_3d, IndexDn_3d)
OppSsqr =NewOperator("Ssqr" ,NF, IndexUp_3d, IndexDn_3d)
OppSplus=NewOperator("Splus",NF, IndexUp_3d, IndexDn_3d)
OppSmin =NewOperator("Smin" ,NF, IndexUp_3d, IndexDn_3d)

OppLx   =NewOperator("Lx"   ,NF, IndexUp_3d, IndexDn_3d)
OppLy   =NewOperator("Ly"   ,NF, IndexUp_3d, IndexDn_3d)
OppLz   =NewOperator("Lz"   ,NF, IndexUp_3d, IndexDn_3d)
OppLsqr =NewOperator("Lsqr" ,NF, IndexUp_3d, IndexDn_3d)
OppLplus=NewOperator("Lplus",NF, IndexUp_3d, IndexDn_3d)
OppLmin =NewOperator("Lmin" ,NF, IndexUp_3d, IndexDn_3d)

OppJx   =NewOperator("Jx"   ,NF, IndexUp_3d, IndexDn_3d)
OppJy   =NewOperator("Jy"   ,NF, IndexUp_3d, IndexDn_3d)
OppJz   =NewOperator("Jz"   ,NF, IndexUp_3d, IndexDn_3d)
OppJsqr =NewOperator("Jsqr" ,NF, IndexUp_3d, IndexDn_3d)
OppJplus=NewOperator("Jplus",NF, IndexUp_3d, IndexDn_3d)
OppJmin =NewOperator("Jmin" ,NF, IndexUp_3d, IndexDn_3d)

Oppldots=NewOperator("ldots",NF, IndexUp_3d, IndexDn_3d)

-- as in the previous example we define the Coulomb interaction

OppF0 =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {1,0,0})
OppF0.Name="F[0]"
OppF2 =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,1,0})
OppF2.Name="F[2]"
OppF4 =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,0,1})
OppF4.Name="F[4]"

-- as in the previous example we define the crystal-field operator

function Akmxyz(Eeg, Et2g)
return ({{0, 0, (2/5)*(Eeg) + (3/5)*(Et2g)} , 
         {4, 0, (21/10)*(Eeg + (-1)*(Et2g))} , 
         {4,-4, (3/2)*((sqrt(7/10))*(Eeg + (-1)*(Et2g)))} , 
         {4, 4, (3/2)*((sqrt(7/10))*(Eeg + (-1)*(Et2g)))} })
end

Eeg  =  0.6
Et2g = -0.4
OpptenDq = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akmxyz(Eeg, Et2g))

-- and an operator that can count the eg and t2g occupation

OppNeg  = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akmxyz(1, 0));
OppNeg.Name="N_eg"
OppNt2g = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akmxyz(0, 1));
OppNt2g.Name="N_t2g"

-- and the core - valence coulomb interaction

OppUsdF0 = NewOperator("U", NF, IndexUp_1s, IndexDn_1s, IndexUp_3d, IndexDn_3d, {1}, {0})
OppUsdG2 = NewOperator("U", NF, IndexUp_1s, IndexDn_1s, IndexUp_3d, IndexDn_3d, {0}, {1})

OppUpdF0 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {1,0}, {0,0})
OppUpdF2 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,1}, {0,0})
OppUpdG1 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {1,0})
OppUpdG3 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {0,1})

Opp2pldots= NewOperator("ldots", NF, IndexUp_2p, IndexDn_2p)

-- and the dipole transition matrices between the 1s and 2p orbital

Akm = {{1,-1,sqrt(1/2)},{1, 1,-sqrt(1/2)}}
TXESx = NewOperator("CF", NF, IndexUp_1s, IndexDn_1s, IndexUp_2p, IndexDn_2p, Akm)
-- y polarized light is defined as y = Sin[phi]Sin[theta] = sqrt(1/2) I ( C_1^{(-1)} + C_1^{(1)})
Akm = {{1,-1,sqrt(1/2)*I},{1, 1,sqrt(1/2)*I}}
TXESy = NewOperator("CF", NF, IndexUp_1s, IndexDn_1s, IndexUp_2p, IndexDn_2p, Akm)
-- z polarized light is defined as z = Cos[theta] = C_1^{(0)}
Akm = {{1,0,1}}
TXESz = NewOperator("CF", NF, IndexUp_1s, IndexDn_1s, IndexUp_2p, IndexDn_2p, Akm)

-- and the annihilation operator on the 1s orbital.

An1sUp = NewOperator("An", NF, IndexUp_1s,{1})
An1sDn = NewOperator("An", NF, IndexDn_1s,{1})

-- once all operators are defined we can set some parameter values.

-- the value of U drops out of a crystal-field calculation as the total number of electrons
-- is always the same
U       =  0.000 
-- F2 and F4 are often referred to in the literature as J_{Hund}. They represent the energy
-- differences between different multiplets. Numerical values can be found in the back of
-- my PhD. thesis for example. http://arxiv.org/abs/cond-mat/0505214 
F2dd    = 11.142 
F4dd    =  6.874
-- F0 is not the same as U, although they are related. Unimportant in crystal-field theory
-- the difference between U and F0 is so important that I do include it here. (U=0 so F0 is not)
F0dd    = U+(F2dd+F4dd)*2/63
-- The Slater integrals between the 2p and 3d shell, again the numerical values can be found
-- in the back of my PhD. thesis. (http://arxiv.org/abs/cond-mat/0505214)
F2pd    =  6.667
G1pd    =  4.922
G3pd    =  2.796
-- in crystal field theory U drops out of the equation, also true for the interaction between the 
-- Ni 1s and Ni 3d electrons
Usd     =  0.000 
-- The Slater integrals between the 1s and 3d shell, can be calculated from the Radial wave
-- functions. You can find these with Cowans code, FPLO, or the Mathematica version of Quanty
G2sd    =  0.057
-- F0 is not the same as U, although they are related. Unimportant in crystal-field theory
-- the difference between U and F0 is so important that I do include it here. (U=0 so F0 is not)
F0sd    =  Usd + G2sd*1/10
-- tenDq in NiO is 1.1 eV as can be seen in optics or using IXS to measure d-d excitations
tenDq   =  1.100
-- the Ni 3d spin-orbit is small but finite
zeta_3d =  0.081
-- the Ni 2p spin-orbit is very large and should not be scaled as theory is quite accurate here
zeta_2p = 11.498
-- we can add a small magnetic field, just to get nice expectation values. (units in eV... )
Bz      = 0.000001
-- In mean field theory the neighboring Ni sites give an effective potential acting on the
-- spin only when magnetically ordered. This exchange field in NiO is 6 J with J=27 meV.
Hex     = 6*0.027 -- see Europhys. Lett., 32 259 (1995) [ and Phys. Rev. B 82, 094403 (2010) ]

-- the total Hamiltonian is the sum of the different operators multiplied with their prefactor
Hamiltonian = F0dd*OppF0 + F2dd*OppF2 + F4dd*OppF4 + tenDq*OpptenDq + zeta_3d*Oppldots + Bz*(2*OppSz+OppLz) + Hex*sqrt(1/6)*(OppSx+OppSy+2*OppSz)
Hamiltonian.Name="H"

-- We normally do not include core-valence interactions between filed and partially filled 
-- shells as it drops out anyhow. For the XAS we thus need to define a "different" 
-- (more complete) Hamiltonian.
XASHamiltonian = Hamiltonian + F0sd * OppUsdF0 + G2sd * OppUsdG2

XESHamiltonian = Hamiltonian + zeta_2p * Opp2pldots + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3


-- We saw in the previous example that NiO has a ground-state triplet with occupation 
-- t2g^6 eg^2 and S=1 (S^2=S(S+1)=2). The next state is 1.1 eV higher in energy and thus
-- unimportant for the ground-state upto temperatures of 10 000 Kelvin. We thus restrict 
-- the calculation to the lowest 3 eigenstates.
Npsi=3
-- in order to make sure we have a filling of 8
-- electrons we need to define some restrictions
-- We need to restrict the occupation of the Ni-1s shell to 2 for the ground state and for
-- the Ni 3d-shell to 8.
StartRestrictions = {NF, NB, {"11 111111 0000000000",8,8}, {"00 000000 1111111111",8,8}}

-- And calculate the lowest 3 eigenfunctions
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)

-- In order to get some information on these eigenstates it is good to plot expectation values
-- We first define a list of all the operators we would like to calculate the expectation value of
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, Oppldots, OppF2, OppF4, OppNeg, OppNt2g};

-- we print the expectation values
PrintExpectationValues(psiList,oppList,{{"zero",0.001}})

-- we here assume that we can create the 1s core hole and that afterwards this core hole 
-- directly decays via the 2p to 1s FY decay. There is no internal relaxation. 
psiListIntermediate = {An1sDn * psiList[1], An1sUp * psiList[1],
                       An1sDn * psiList[2], An1sUp * psiList[2],
                       An1sDn * psiList[3], An1sUp * psiList[3]}


-- calculating the spectra is simple and single line once all operators and wave-functions
-- are defined.
XESSpectra = CreateSpectra(XESHamiltonian, {TXESx,TXESy,TXESz}, psiListIntermediate, {{"Emin",-10}, {"Emax",25}, {"NE",3500}, {"Gamma",0.2}})
-- XASSpectra will contain 15 spectra in the end. We have 3 functions in psiList and 5 
-- polarizations. We furthermore define the energy range to be from -5 to +10
-- using 3000 energy points (delta E = 10 meV) and an internal lifetime during the calculation
-- of 100 meV (Full width Half maximum)

-- we can use several operations on a spectrum object. For one we can sum several spectra
-- the function Spectra.Sum sums the spectra with the prefactors given in the list afterwards
-- sum the spectra of the ground-state for different polarizations in order to get the
-- isotropic spectrum
XESIsoSpectra = Spectra.Sum(XESSpectra,{1,1,0,0,0,0,  1,1,0,0,0,0, 1,1,0,0,0,0})
-- broaden the spectra, note that a method changes the spectrum object, here we broaden by
-- a Gaussian of FWHM 0.4 and Lorenzian of FWHM 1.0
XESIsoSpectra.Broaden(0.4,1.0)

-- in order to plot the spectra we write them to file in ASCII format
-- note that the calculated object, i.e. <psi | T^dag 1/(w-H+IG/2) T | psi>
-- is complex. Minus the imaginary part is the absorption.
XESSpectra.Print({{"file","XESSpec.dat"}})
XESIsoSpectra.Print({{"file","XESIsoSpec.dat"}})

-- in order to plot we use gnuplot. You can use any program, but gnuplot is nice as it can
-- be called from within Quanty. We first define the gnuplot script as a string inside Quanty
gnuplotInput = [[
set autoscale 
set xtic auto 
set ytic auto 
set style line  1 lt 1 lw 1 lc rgb "#000000"
set style line  2 lt 1 lw 1 lc rgb "#FF0000"
set style line  3 lt 1 lw 3 lc rgb "#000000"

set xlabel "E (eV)" font "Times,12"
set ylabel "Intensity (arb. units)" font "Times,8"

set out 'XESSpec.ps'
set size 1.0, 1.0
set terminal postscript landscape enhanced color  "Times" 8

energyshift=7475

set multiplot layout 6, 3

plot "XESSpec.dat"  u (-$1+energyshift):(-$3 ) title 'x-polarized dn Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$5 ) title 'x-polarized up Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$7 ) title 'x-polarized dn Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$9 ) title 'x-polarized up Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$11) title 'x-polarized dn Sz= 1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$13) title 'x-polarized up Sz= 1' with lines ls  1

plot "XESSpec.dat"  u (-$1+energyshift):(-$15) title 'y-polarized dn Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$17) title 'y-polarized up Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$19) title 'y-polarized dn Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$21) title 'y-polarized up Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$23) title 'y-polarized dn Sz= 1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$25) title 'y-polarized up Sz= 1' with lines ls  1

plot "XESSpec.dat"  u (-$1+energyshift):(-$27) title 'z-polarized dn Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$29) title 'z-polarized up Sz=-1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$31) title 'z-polarized dn Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$33) title 'z-polarized up Sz= 0' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$35) title 'z-polarized dn Sz= 1' with lines ls  1
plot "XESSpec.dat"  u (-$1+energyshift):(-$37) title 'z-polarized up Sz= 1' with lines ls  1


unset multiplot


set size 1.0, 0.6
intensityscale=48
set out 'XESIsoSpec.ps'

plot "XESIsoSpec.dat" using (-$1+energyshift):((-$3) * intensityscale) title 'isotropic theory' with lines ls 3

]]

-- next we save this string to a file
file = io.open("XESSpec.gnuplot", "w")
file:write(gnuplotInput)
file:close()

-- and finally call gnuplot to execute the script
os.execute("gnuplot XESSpec.gnuplot")
-- as I like pdf to view I transform the format
os.execute(" ps2pdf XESSpec.ps ; ps2pdf XESIsoSpec.ps")