-- This tutorial calculates the 3d inverse photo emission spectra of Ni in NiO using
-- crystal field theory

-- An electron from the a free electron is excited into the Ni 3d shell producing a photon. We can think
-- of these incoming electrons 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, Index, NF, epsilon_x, epsilon_y, epsilon_z, detector_theta, detector_phi)
  local T = 0
  local epsilon={}
  local l1 = (#Index - 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("Cr", NF, {Index[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

-- 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. In order to calculate 2p to 3d excitations we
-- also need a Ni 2p shell. We thus have a total of 10+6=16 fermions, 6 Ni-2p and 10 Ni-3d
-- spin-orbitals
NF=16
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}

-- 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})
OppF2 =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,1,0})
OppF4 =NewOperator("U", NF, IndexUp_3d, IndexDn_3d, {0,0,1})

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

Akm = PotentialExpandedOnClm("Oh",2,{0.6,-0.4})
OpptenDq = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)

-- and as in the previous example we define operators that count the number of eg and t2g
-- electrons

Akm = PotentialExpandedOnClm("Oh",2,{1,0})
OppNeg = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)
Akm = PotentialExpandedOnClm("Oh",2,{0,1})
OppNt2g = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm)
OppN = OppNeg+OppNt2g

-- 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. However the moment we change the number of electrons in PES and IPES
-- we need to set U and ed to get the correct gap
nd      =  8
U       =  7.3 
-- ed is set such that the total energy E = U nd(nd-1)/2 + ed nd has a minimum for nd=8
ed      =  -U*(nd-1/2)
-- 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
-- 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
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) + ed*OppN

-- 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-2p shell to 6 for the ground state and for
-- the Ni 3d-shell to 8.
StartRestrictions = {NF, NB, {"111111 0000000000",6,6}, {"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};

-- next we loop over all operators and all states and print the expectation value
print(" <E>    <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <l.s>  <F[2]> <F[4]> <Neg>  <Nt2g>");
for i = 1,#psiList do
  for j = 1,#oppList do
    expectationvalue = Chop(psiList[i]*oppList[j]*psiList[i])
    io.write(string.format("%6.3f ",Complex.Re(expectationvalue)))
  end
  io.write("\n")
end

epsilon = {1,0,0}
detector_theta = pi/4
detector_phi = 0
TIPES_p = TransitionOperator (1, IndexDn_3d, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TIPES_f = TransitionOperator (3, IndexDn_3d, NF, epsilon[1], epsilon[2], epsilon[3], detector_theta, detector_phi)
TIPES = TIPES_p + TIPES_f

-- calculating the spectra is simple and single line once all operators and wave-functions
-- are defined.
IPESSpectra = CreateSpectra(Hamiltonian, {TIPES_p,TIPES_f,TIPES}, psiList, {{"Emin",-10}, {"Emax",15}, {"NE",2500}, {"Gamma",0.1}})
-- 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 -10 to +15
-- using 2500 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
IPESSpectra.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 inverse photo-emission intensity.
IPESSpectra.Print({{"file","IPESSpectra.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 'IPESSpectra.ps'
set size 1.0, 1.0
set terminal postscript landscape enhanced color  "Times" 12

set multiplot layout 3, 3

plot "IPESSpectra.dat"  u 1:(-$3 ) title '  s Sz=-1' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$5 ) title '  s Sz= 0' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$7 ) title '  s Sz= 1' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$9 ) title '  d Sz=-1' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$11) title '  d Sz= 0' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$13) title '  d Sz= 1' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$15) title 's+d Sz=-1' with lines ls  1
plot "IPESSpectra.dat"  u 1:(-$17) title 's+d Sz= 0' with lines ls  1
plot "IPESSpectra.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("IPESSpectra.gnuplot", "w")
file:write(gnuplotInput)
file:close()

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