-- This tutorial calculates the 3d x-ray photo emission spectra of Ni in NiO using
-- an Anderson impurity model (occupied bath sites only).  We build upon the ligand field
-- theory example, and use the function in
dofile("../Expand_Ligand_Basis.Quanty")
-- to expand the basis


-- We define the ligand field theory basis

NF=20
NB=0
IndexDn_3d={ 0, 2, 4, 6, 8}
IndexUp_3d={ 1, 3, 5, 7, 9}
IndexDn_Ld={10,12,14,16,18}
IndexUp_Ld={11,13,15,17,19}

-- 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)
U       =  7.3
Delta   =  4.7
-- parameters obtained from DFT (PRB 85, 165113 (2012))
F2dd    = 11.142 
F4dd    =  6.874
tenDq   =  0.56
tenDqL  =  1.44
Veg     =  2.06
Vt2g    =  1.21
zeta_3d =  0.081
Bz      =  0.000001

-- turning U and Delta to onsite energies (Including the transformation from U to F0)
-- Also including the shift of zero energy to account for the chemical potential

ed      = (10*Delta-nd*(19+nd)*U/2)/(10+nd)
eL      = nd*((1+nd)*U/2-Delta)/(10+nd)
e       = - ed - (2*nd - 1)*U/2
eshift  = -3

edPrime = ed + e + eshift
eLPrime = eL + e + eshift

F0dd    = U+(F2dd+F4dd)*2/63

-- 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
eLPrime      = eLPrime + eLshift
epsilon0Bath = eLPrime * Matrix.Identity(5)
tMat         = (W/4)   * Matrix.Identity(5)

EshiftTotal = 10*NShells*eLPrime

-- Define bit masks that can be used to select different shells, on the original basis
BitMasks = {"1111111111 0000000000", "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],nd,nd}, {BitMasks[2],10,10}, {BitMasks[3],10*NShells,10*NShells}}
GroundStateRestrictions = {NF, NB, {BitMasks[1],nd,nd+dExtraMax}, {BitMasks[3],10*NShells-HoleBathMax,10*NShells}}
PESRestrictions         = {NF, NB, {BitMasks[1],nd-1,nd+dExtraMax}, {BitMasks[3],10*NShells-HoleBathMax,10*NShells}}

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



-- An electron from the 3d shell is excited by a photon to a free electron. We can think
-- of these as plane-waves, but equally well as spherical harmonics times Bessel functions
-- For the later it becomes clear that the angular momentum of the free electron is either
-- 1 (p orbital) or 3 (f orbital). The ratio between these two channels depends 
-- on the photon energy
-- There is furthermore an angular dependence. i.e. depending on the polarisation and 
-- direction of the detector one finds a different intensity  
-- Here we define the transition operator, depending on the polarisation and detector position

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={}
  local l1 = (#IndexDn - 1)/2
  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=-l1,l1 do -- loop over angular momentum of the bound 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+l1+1]}, {epsilon[m2] * YlmClmYlm(l1,m1,1,m2,l3,m3) * SphericalHarmonicY(l3,m3,detector_theta,detector_phi)})
              + NewOperator("An", NF, {IndexUp[m1+l1+1]}, {epsilon[m2] * YlmClmYlm(l1,m1,1,m2,l3,m3) * SphericalHarmonicY(l3,m3,detector_theta,detector_phi)})
      end
    end
  end
  return Chop(T)
end


-- we can define the angular momentum operators for the d-shell as we did in crystal-field
-- theory

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)

-- And similar we can define the angular momentum operators for the ligand d-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)

-- In order to calculate the total angular momentum in the cluster we can sum the operators

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

-- Just like in crystal-field theory we have Coulomb interaction on the Ni-d shell.
-- We again expand the Coulomb interaction on spherical harmonics, where the angular part
-- is solved analytical and the radial part gives three parameters F0, F2 and F4.
-- We here define three operators separately and only later provide parameters

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})

-- In ligand-field theory the ligand-field interaction is given by three different terms.
-- There is an onsite splitting on the Transition metal d-shell
-- There is an onsite splitting on the Ligand d-shell
-- There is a hopping between the ligand d-shell and the Transition metal d-shell

-- These interactions can be seen as effective potentials responsible for the splitting
-- In order to enter these potentials we expand them on renormalized spherical harmonics
-- and add the expansion coefficients to the function NewOperator("CF", ...)
-- We thus need to know the potential expanded on spherical harmonics: 
-- Akm = {{k1,m1,Akm1},{k2,m2,Akm2}, ... }

-- For specific symmetries we can use the function "PotentialExpandedOnClm" Which for cubic
-- symmetry needs the energy of the eg and t2g orbitals. We here take the potential to be
-- such that we have a 1 eV splitting and later multiply the operator with the actual size

-- In crystal-field theory there is only an interaction on the transition metal d-shell
-- In ligand field theory there is an interaction on the transition metal d-shell as well
-- as on the ligand d-shell

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)

-- We want to be able to calculate the occupation of the eg and t2g orbitals, we here use
-- the same operators with potentials of 1 for the eg or 1 for the t2g orbitals to create
-- number operators. (Note that there are many other options to do this)

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)

-- We also want to know hom many electrons are in the Ni-d and how many are in the Ligand-d 
-- shell. Here the number operators that count them.

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

-- Besides the onsite energy of the ligand and transition metal d-shell we need to define the 
-- hopping between them. We can use the same crystal-field operator, but now acting between
-- two different shells.

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)








-- 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)
            + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g
            + edPrime * OppN_3d + eLPrime * OppN_Ld
            + 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

epsilon = {1,0,0}
detector_theta = pi/4
detector_phi = 0
TPES_p = TransitionOperator (1, IndexDn_3d, IndexUp_3d, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TPES_f = TransitionOperator (3, IndexDn_3d, IndexUp_3d, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TPES = TPES_p + TPES_f

-- calculating the spectra is simple and single line once all operators and wave-functions
-- are defined.
PESSpectra = CreateSpectra(Hamiltonian, {TPES_p,0*TPES_f,0*TPES}, psiList, {{"Emin",-5}, {"Emax",25}, {"NE",3000}, {"Gamma",0.1}, {"Restrictions",PESRestrictions}})
-- PESSpectra will contain 9 spectra in the end. We have 3 functions in psiList and 3 
-- transition operators. We furthermore define the energy range to be from -5 to +25
-- using 3000 energy points (delta E = 10 meV) and an internal lifetime during the calculation
-- of 100 meV (Full width Half maximum)

-- 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
PESSpectra.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 photo-emission intensity.
PESSpectra.Print({{"file","PESSpectra.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,12"

set out 'PESSpectra.ps'
set size 1.0, 1.0
set terminal postscript landscape enhanced color  "Times" 12

set multiplot layout 3, 3

plot "PESSpectra.dat"  u 1:(-$3 ) title '  s Sz=-1' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$5 ) title '  s Sz= 0' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$7 ) title '  s Sz= 1' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$9 ) title '  d Sz=-1' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$11) title '  d Sz= 0' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$13) title '  d Sz= 1' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$15) title 's+d Sz=-1' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$17) title 's+d Sz= 0' with lines ls  1
plot "PESSpectra.dat"  u 1:(-$19) title 's+d Sz= 1' with lines ls  1

unset multiplot

]]

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

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