-- this example calculates the resonant inelastic x-ray scattering in NiO we look at the
-- Ni L23M45 edge, i.e. we make an excitation from 2p to 3d (L23) and decay from the
-- 3d shell back to the 2p shell (final "hole" in the 3d shell M45 edge). These spectra
-- measure d-d excitations and or magnons.

-- The calculations are done using an Anderson impurity model, whereby we start from a 
-- Ligand field theory input and make minor modifications
dofile("../Expand_Ligand_Basis.Quanty")

Verbosity(0)
-- This tutorial calculates the L23 XSAS of NiO within the Anderson impurity approximation.
-- We start with a ligand field theory model and then replace the single bath site by a band
-- Have a look at the ligand field theory calculations first.

-- We define the ligand field theory basis
 
NF=26
NB=0
IndexDn_2p={ 0, 2, 4}
IndexUp_2p={ 1, 3, 5}
IndexDn_3d={ 6, 8,10,12,14}
IndexUp_3d={ 7, 9,11,13,15}
IndexDn_Ld={16,18,20,22,24}
IndexUp_Ld={17,19,21,23,25}
 
-- We bring the parameters to the front in order to be able to extend the bath

-- number of electrons (formal valence)
nd = 8
-- parameters from experiment (core level PES)
Udd     =  7.3
Upd     =  8.5
Delta   =  4.7
-- parameters obtained from DFT (PRB 85, 165113 (2012))
F2dd    = 11.14 
F4dd    =  6.87
F2pd    =  6.67
G1pd    =  4.92
G3pd    =  2.80
tenDq   =  0.56
tenDqL  =  1.44
Veg     =  2.06
Vt2g    =  1.21
zeta_3d =  0.081
zeta_2p = 11.51
Bz      =  0.000001
Hz      =  0.120
 
ed      = (10*Delta-nd*(19+nd)*Udd/2)/(10+nd)
eL      = nd*((1+nd)*Udd/2-Delta)/(10+nd)
 
epfinal = (10*Delta + (1+nd)*(nd*Udd/2-(10+nd)*Upd)) / (16+nd)
edfinal = (10*Delta - nd*(31+nd)*Udd/2-90*Upd) / (16+nd)
eLfinal = ((1+nd)*(nd*Udd/2+6*Upd) - (6+nd)*Delta) / (16+nd)
 
F0dd    = Udd + (F2dd+F4dd) * 2/63
F0pd    = Upd + (1/15)*G1pd + (3/70)*G3pd

-- additional parameters for the bath.  We assume a semicircular O-2p DOS. One can make this
-- much more realistic, using the DFT bath, More on this in later tutorials.

-- the band width
W = 6
-- the number of discretisation points
NShells = 10 -- this yields 6 ligand shells
-- restrictions
HoleBathMax   = 1
dExtraMax     = 2
-- additional shift of all bath energies, including the ligand shell needed while
-- although the hopping goes to the bath center, the bath edge have stronger hybridization
eLshift = -1.5

-- derived parameters
eL           = eL + eLshift
epsilon0Bath = eL * Matrix.Identity(5)
tMat         = (W/4) * Matrix.Identity(5)

EshiftTotal = 10*NShells*eL

-- Define bit masks that can be used to select different shells, on the original basis
BitMasks = {"111111 0000000000 0000000000", "000000 1111111111 0000000000", "000000 0000000000 1111111111"}
-- Expand the basis, using the function defined in the file ../ExpandLigandBasis.Quanty
NF, HBath, BitMasks = ExpandLigandBasis(NShells, NF, IndexDn_Ld, IndexUp_Ld, tMat, epsilon0Bath, BitMasks)


-- Define restrictions
StartRestrictions       = {NF, NB, {BitMasks[1],6,6}, {BitMasks[2],nd,nd}, {BitMasks[3],10,10}, {BitMasks[4],10*NShells,10*NShells}}
GroundStateRestrictions = {NF, NB, {BitMasks[2],nd,nd+dExtraMax}, {BitMasks[4],10*NShells-HoleBathMax,10*NShells}}
XASRestrictions         = GroundStateRestrictions



-- the script now follows the ligand field script for a large part

-- angular momentum operators on the d-shell
 
OppSx_3d   =NewOperator("Sx"   ,NF, IndexUp_3d, IndexDn_3d)
OppSy_3d   =NewOperator("Sy"   ,NF, IndexUp_3d, IndexDn_3d)
OppSz_3d   =NewOperator("Sz"   ,NF, IndexUp_3d, IndexDn_3d)
OppSsqr_3d =NewOperator("Ssqr" ,NF, IndexUp_3d, IndexDn_3d)
OppSplus_3d=NewOperator("Splus",NF, IndexUp_3d, IndexDn_3d)
OppSmin_3d =NewOperator("Smin" ,NF, IndexUp_3d, IndexDn_3d)
 
OppLx_3d   =NewOperator("Lx"   ,NF, IndexUp_3d, IndexDn_3d)
OppLy_3d   =NewOperator("Ly"   ,NF, IndexUp_3d, IndexDn_3d)
OppLz_3d   =NewOperator("Lz"   ,NF, IndexUp_3d, IndexDn_3d)
OppLsqr_3d =NewOperator("Lsqr" ,NF, IndexUp_3d, IndexDn_3d)
OppLplus_3d=NewOperator("Lplus",NF, IndexUp_3d, IndexDn_3d)
OppLmin_3d =NewOperator("Lmin" ,NF, IndexUp_3d, IndexDn_3d)
 
OppJx_3d   =NewOperator("Jx"   ,NF, IndexUp_3d, IndexDn_3d)
OppJy_3d   =NewOperator("Jy"   ,NF, IndexUp_3d, IndexDn_3d)
OppJz_3d   =NewOperator("Jz"   ,NF, IndexUp_3d, IndexDn_3d)
OppJsqr_3d =NewOperator("Jsqr" ,NF, IndexUp_3d, IndexDn_3d)
OppJplus_3d=NewOperator("Jplus",NF, IndexUp_3d, IndexDn_3d)
OppJmin_3d =NewOperator("Jmin" ,NF, IndexUp_3d, IndexDn_3d)
 
Oppldots_3d=NewOperator("ldots",NF, IndexUp_3d, IndexDn_3d)
 
-- Angular momentum operators on the Ligand shell
 
OppSx_Ld   =NewOperator("Sx"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppSy_Ld   =NewOperator("Sy"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppSz_Ld   =NewOperator("Sz"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppSsqr_Ld =NewOperator("Ssqr" ,NF, IndexUp_Ld, IndexDn_Ld)
OppSplus_Ld=NewOperator("Splus",NF, IndexUp_Ld, IndexDn_Ld)
OppSmin_Ld =NewOperator("Smin" ,NF, IndexUp_Ld, IndexDn_Ld)
 
OppLx_Ld   =NewOperator("Lx"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppLy_Ld   =NewOperator("Ly"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppLz_Ld   =NewOperator("Lz"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppLsqr_Ld =NewOperator("Lsqr" ,NF, IndexUp_Ld, IndexDn_Ld)
OppLplus_Ld=NewOperator("Lplus",NF, IndexUp_Ld, IndexDn_Ld)
OppLmin_Ld =NewOperator("Lmin" ,NF, IndexUp_Ld, IndexDn_Ld)
 
OppJx_Ld   =NewOperator("Jx"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppJy_Ld   =NewOperator("Jy"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppJz_Ld   =NewOperator("Jz"   ,NF, IndexUp_Ld, IndexDn_Ld)
OppJsqr_Ld =NewOperator("Jsqr" ,NF, IndexUp_Ld, IndexDn_Ld)
OppJplus_Ld=NewOperator("Jplus",NF, IndexUp_Ld, IndexDn_Ld)
OppJmin_Ld =NewOperator("Jmin" ,NF, IndexUp_Ld, IndexDn_Ld)
 
-- total angular momentum
OppSx = OppSx_3d + OppSx_Ld
OppSy = OppSy_3d + OppSy_Ld
OppSz = OppSz_3d + OppSz_Ld
OppSsqr = OppSx * OppSx + OppSy * OppSy + OppSz * OppSz
OppLx = OppLx_3d + OppLx_Ld
OppLy = OppLy_3d + OppLy_Ld
OppLz = OppLz_3d + OppLz_Ld
OppLsqr = OppLx * OppLx + OppLy * OppLy + OppLz * OppLz
OppJx = OppJx_3d + OppJx_Ld
OppJy = OppJy_3d + OppJy_Ld
OppJz = OppJz_3d + OppJz_Ld
OppJsqr = OppJx * OppJx + OppJy * OppJy + OppJz * OppJz
 
-- define the coulomb operator
 
OppF0_3d =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {1,0,0})
OppF2_3d =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,1,0})
OppF4_3d =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,0,1})
 
-- crystal and ligand field
 
Akm = PotentialExpandedOnClm("Oh", 2, {0.6,-0.4})
OpptenDq_3d = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)
OpptenDq_Ld = NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, Akm)
 
Akm = PotentialExpandedOnClm("Oh", 2, {1,0})
OppNeg_3d = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)
OppNeg_Ld = NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, Akm)
Akm = PotentialExpandedOnClm("Oh", 2, {0,1})
OppNt2g_3d = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)
OppNt2g_Ld = NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, Akm)
 
OppNUp_2p = NewOperator("Number", NF, IndexUp_2p, IndexUp_2p, {1,1,1})
OppNDn_2p = NewOperator("Number", NF, IndexDn_2p, IndexDn_2p, {1,1,1})
OppN_2p = OppNUp_2p + OppNDn_2p
OppNUp_3d = NewOperator("Number", NF, IndexUp_3d, IndexUp_3d, {1,1,1,1,1})
OppNDn_3d = NewOperator("Number", NF, IndexDn_3d, IndexDn_3d, {1,1,1,1,1})
OppN_3d = OppNUp_3d + OppNDn_3d
OppNUp_Ld = NewOperator("Number", NF, IndexUp_Ld, IndexUp_Ld, {1,1,1,1,1})
OppNDn_Ld = NewOperator("Number", NF, IndexDn_Ld, IndexDn_Ld, {1,1,1,1,1})
OppN_Ld = OppNUp_Ld + OppNDn_Ld
 
-- define L-d interaction
 
Akm = PotentialExpandedOnClm("Oh", 2, {1,0})
OppVeg  = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) +  NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm)
Akm = PotentialExpandedOnClm("Oh", 2, {0,1})
OppVt2g = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) +  NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm)
 
-- core valence interaction
 
Oppcldots= NewOperator("ldots", NF, IndexUp_2p, IndexDn_2p)
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})
 
-- dipole transition
 
t=math.sqrt(1/2)
 
Akm = {{1,-1,t},{1, 1,-t}}
T2p3dx = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
Akm = {{1,-1,t*I},{1, 1,t*I}}
T2p3dy = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
Akm = {{1,0,1}}
T2p3dz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 
T2p3dr = t*(T2p3dx - I * T2p3dy)
T2p3dl =-t*(T2p3dx + I * T2p3dy)

-- we can remove zero's from the dipole operator by chopping it.
T2p3dr.Chop()
T2p3dl.Chop()

-- the 3d to 2p dipole transition is the conjugate transpose of the 2p to 3d dipole transition
T3d2px = ConjugateTranspose(T2p3dx)
T3d2py = ConjugateTranspose(T2p3dy)
T3d2pz = ConjugateTranspose(T2p3dz)
T3d2pl = ConjugateTranspose(T2p3dl)
T3d2pr = ConjugateTranspose(T2p3dr)
 
-- and our Hamiltonian is the sum over several operators
 
Hamiltonian    =  F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d) + Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + ed * OppN_3d + eL * OppN_Ld
               + HBath - EshiftTotal
XASHamiltonian =  F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d)+ Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + edfinal * OppN_3d + eLfinal * OppN_Ld + epfinal * OppN_2p + zeta_2p * Oppcldots + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3  
               + HBath - EshiftTotal
            
-- we now can create the lowest Npsi eigenstates:
Npsi=3

psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi,{{"Restrictions",GroundStateRestrictions}})
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz_3d, OppLz_3d, Oppldots_3d, OppF2_3d, OppF4_3d, OppNeg_3d, OppNt2g_3d, OppNeg_Ld, OppNt2g_Ld, OppN_3d}
 
-- print of some expectation values
 
print("  #    <E>      <S^2>    <L^2>    <J^2>    <S_z^3d> <L_z^3d> <l.s>    <F[2]>   <F[4]>   <Neg^3d> <Nt2g^3d><Neg^Ld> <Nt2g^Ld><N^3d>");
for i = 1,#psiList do
  io.write(string.format("%3i ",i))
  for j = 1,#oppList do
    expectationvalue = Chop(psiList[i]*oppList[j]*psiList[i])
    io.write(string.format("%8.3f ",expectationvalue))
  end
  io.write("\n")
end


-- Not the main task of this example, but in order to know where the resonance sits it is
-- nice to calculate the x-ray absorption spectra
XASSpectra = CreateSpectra(XASHamiltonian, T2p3dx, psiList[1], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0},{"Restrictions",XASRestrictions}})
XASSpectra.Print({{"file", "RIXSL23M45_XAS.dat"}})

-- We also calculate the fluorescence yield spectra which are equal to the integral over the
-- RIXS spectra
FYSpectra = CreateFluorescenceYield(XASHamiltonian, T2p3dx, T3d2py, psiList[1], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0},{"Restrictions",XASRestrictions}})
FYSpectra.Print({{"file","RIXSL23M45_FY.dat"}})

-- and we calculate the RIXS spectra. Note that in order to calculate RIXS you need to
-- specify two Hamiltonians.
-- These calculations can be lengthy and one should think that for each incoming energy
-- (E1) one needs to do a new calculation. In this case there are thus 120 calculations
RIXSSpectra = CreateResonantSpectra(XASHamiltonian, Hamiltonian, T2p3dx, T3d2py, psiList[1], {{"Emin1",-10}, {"Emax1",20}, {"NE1",120}, {"Gamma1",1.0}, {"Emin2",-0.5}, {"Emax2",7.5}, {"NE2",8000}, {"Gamma2",0.05},{"Restrictions1",XASRestrictions},{"Restrictions2",GroundStateRestrictions},{"DenseBorder",10000}})
RIXSSpectra.Print({{"file", "RIXSL23M45.dat"}})

print("Finished calculating the spectra now start plotting.\nThis might take more time than the calculation");

-- Once finished we can make some nice plots. The spectra are saved to disk in ASCII format
-- but I like to add gnuplot scripts so you can look at pictures imeidately
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 1 lc rgb "#00FF00"
set style line  4 lt 1 lw 1 lc rgb "#0000FF"

set out 'RIXSL23M45_Map.ps'
set terminal postscript portrait enhanced color  "Times" 8 size 7.5,6

unset colorbox

energyshift=857.6

set multiplot 
set size 0.5,0.55
set origin 0,0

set ylabel "resonant energy (eV)" font "Times,10"
set xlabel "energy loss (eV)" font "Times,10"

set yrange [852:860]
set xrange [-0.5:7.5]

plot "<(awk '((NR>5)&&(NR<8007)){for(i=3;i<=NF;i=i+2){printf \"%s \", $i}{printf \"%s\", RS}}' RIXSL23M45.dat)" matrix using ($2/1000-0.5):($1/4+energyshift-10):(-$3) with image notitle 

set origin 0.5,0

set yrange [869:877]
set xrange [-0.5:7.5]

plot "<(awk '((NR>5)&&(NR<8007)){for(i=3;i<=NF;i=i+2){printf \"%s \", $i}{printf \"%s\", RS}}' RIXSL23M45.dat)" matrix using ($2/1000-0.5):($1/4+energyshift-10):(-$3) with image notitle 

unset multiplot

set out 'RIXSL23M45_Spec.ps'
set terminal postscript portrait enhanced color  "Times" 8 size 7.5,5

set multiplot 
set size 0.25,1.0
set origin 0,0

set ylabel "E (eV)" font "Times,10"
set xlabel "XAS" font "Times,10"
set yrange [energyshift-10:energyshift+20]
set xrange [-0.3:0]
plot "RIXSL23M45_XAS.dat"  using       3:($1+energyshift) notitle with lines ls  1,\
     "RIXSL23M45_FY.dat"   using (-3*$2):($1+energyshift) notitle with lines ls  4

set size 0.8,1.0
set origin 0.2,0.0

set xlabel "Energy loss (eV)" font "Times,10"
unset ylabel 
unset ytics
set xrange [-0.5:7.5]

ofset = 0.25 
scale=3

plot for [i=0:120] "RIXSL23M45.dat"  using 1:(-column(3+2*i)*scale+ofset*i-10 + energyshift) notitle with lines ls  4

unset multiplot
]]

-- write the gnuplot script to a file
file = io.open("RIXSL23M45.gnuplot", "w")
file:write(gnuplotInput)
file:close()

-- call gnuplot to execute the script
os.execute("gnuplot RIXSL23M45.gnuplot")
-- transform to pdf and eps
os.execute("ps2pdf RIXSL23M45_Map.ps")
os.execute("ps2pdf RIXSL23M45_Spec.ps")




