
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 = 5 -- 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}}
TXASx = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
Akm = {{1,-1,t*I},{1, 1,t*I}}
TXASy = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
Akm = {{1,0,1}}
TXASz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 
TXASr = t*(TXASx - I * TXASy)
TXASl =-t*(TXASx + I * TXASy)
 
-- 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
 
-- We calculate the x-ray absorption spectra for z, right circular and left circular polarized light for the 3 lowest eigen-states. (9 spectra in total)
 
XASSpectra = CreateSpectra(XASHamiltonian, {TXASz, TXASr, TXASl}, psiList, {{"Emin",-15}, {"Emax",25}, {"NE",2000}, {"Gamma",0.1}, {"Restrictions",XASRestrictions}})
 
-- and put some additional energy broadening on it (0.4 Gaussian and energy dependent lorenzian)
XASSpectra.Broaden(0.4, {{-3.7, 0.45}, {-2.2, 0.65}, { 0.0, 0.65}, { 1.0, 2.00}, { 6  , 2.00}, { 8  , 0.80}, {13.2, 0.80}, {14.0, 0.90}, {16.0, 0.90}, {17.0, 2.00}})
 
-- 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
XASIsoSpectra = Spectra.Sum(XASSpectra,{1,0,0, 1,0,0, 1,0,0})
 
 
-- We can print the spectra to file
XASSpectra.Print({{"file","XASSpec.dat"}})
XASIsoSpectra.Print({{"file","XASIsoSpec.dat"}})
 
-- In order to  make plots I use Gnu-plot, which can be called from within Quanty (if it is installed on your computer)
-- Here a script to make the plots
 
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,12"
 
set out 'XASSpec.ps'
set size 1.0, 1.0
set terminal postscript portrait enhanced color  "Times" 12
 
set multiplot layout 3, 3
 
plot "XASSpec.dat"  u 1:(-$3 ) title 'z-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u 1:(-$5 ) title 'z-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u 1:(-$7 ) title 'z-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u 1:(-$9 ) title 'r-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u 1:(-$11) title 'r-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u 1:(-$13) title 'r-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u 1:(-$15) title 'l-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u 1:(-$17) title 'l-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u 1:(-$19) title 'l-polarized Sz= 1' with lines ls  1
 
unset multiplot
 
energyshift=857.7
intensityscale=64
 
plot "XASSpec.dat" using ($1+energyshift):((-$3-$9-$15) * intensityscale) title 'isotropic theory' with lines ls 1,\
     "../../NiO Experiment/XAS_L23_PRB_57_11623_1998" using 1:2 title 'isotropic experiment' with lines ls 2
 
 
set size 1.0, 0.6
intensityscale=48
set out 'XASIsoSpec.ps'
set xrange [847:877]
 
plot "../../NiO Experiment/XAS_L23_PRB_57_11623_1998" using 1:2 title 'isotropic experiment' with filledcurves y1=0,\
     "XASIsoSpec.dat" using ($1+energyshift):((-$3) * intensityscale) title 'isotropic theory' with lines ls 3
 
]]
 
-- write the gnuplot script to a file
file = io.open("XASSpec.gnuplot", "w")
file:write(gnuplotInput)
file:close()
 
-- and finally call gnuplot to execute the script
os.execute("gnuplot XASSpec.gnuplot")
-- as I like pdf to view and eps to include in the manuel I transform the format
os.execute(" ps2pdf XASSpec.ps ; ps2pdf XASIsoSpec.ps")