-- This tutorial calculates the 1s to 3d x-ray absorption spectra of Ni in NiO using
-- crystal field theory

-- 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=12
NB=0
IndexDn_1s={0}
IndexUp_1s={1}
IndexDn_3d={2,4,6,8,10}
IndexUp_3d={3,5,7,9,11}

-- 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"

-- we also need to define the Coulomb interaction between the Ni 1s- and Ni 3d-shell
-- Again the interaction (e^2/(|r_i-r_j|)) is expanded on spherical harmonics. For the interaction
-- between two shells we need to consider two cases. For the direct interaction a 1s electron
-- scatters of a 3d electron into a 1s and 3d electron. The radial integrals involve
-- the square of a 1s radial wave function at coordinate 1 and the square of a 3d radial
-- wave function at coordinate 2. The transfer of angular momentum must be 0.
-- These processes are called direct and the resulting Slater integrals are F[0].
-- The second proces involves a 1s electron scattering of a 3d electron into the 3d shell
-- and at the same time the 3d electron scattering into a 1s shell. These exchange processes
-- involve radial integrals over the product of a 1s and 3d radial wave function. The transfer
-- of angular momentum in this case must be 2 and the Slater integrals are called G[2].

-- In Quanty you can enter these processes by labeling 4 indices for the orbitals, once
-- the 2p shell with spin up, 2p shell with spin down, 3d shell with spin up and 3d shell with
-- spin down. Followed by the direct Slater integrals (F0) and the exchange Slater 
-- integrals (G2)

-- Here we define the operators separately and later sum them with appropriate prefactors

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

-- next we define the Quadrupole operator. The quadrupole operator is given as epsilon.r * k.r
-- with epsilon the polarization vector of the light, r the unit position vector and
-- k the momentum vector of the light
-- We can expand the position vector on (renormalized) spherical harmonics and use
-- the crystal-field operator to create the quadrupole operator. 
-- Note that a polarisation in the x direction and momentum in the y direction gives the 
-- same operator as for the case where the polarisation is in the y direction and the 
-- momentum is in the x direction

-- There are 5 possible quadrupole operators and all other polarisations and k directions
-- can be made from linear combinations
-- yz "polarisation"
Akm = {{2,-1,I*sqrt(1/2)},{2, 1,I*sqrt(1/2)}}
TXASyz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_1s, IndexDn_1s, Akm)
-- xz "polarisation"
Akm = {{2,-1,  sqrt(1/2)},{2, 1, -sqrt(1/2)}}
TXASxz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_1s, IndexDn_1s, Akm)
-- xy "polarisation"
Akm = {{2,-2,I*sqrt(1/2)},{2, 2,-I*sqrt(1/2)}}
TXASxy = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_1s, IndexDn_1s, Akm)
-- 3z-r^2 "polarisation"
Akm = {{2,0,1}}
TXASz2 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_1s, IndexDn_1s, Akm)
-- x^2-y^2 "polarisation"
Akm = {{2,-2,  sqrt(1/2)},{2, 2,  sqrt(1/2)}}
TXASx2y2 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_1s, IndexDn_1s, Akm)

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

-- 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 0000000000",2,2}, {"00 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}})

-- calculating the spectra is simple and single line once all operators and wave-functions
-- are defined.
XASSpectra = CreateSpectra(XASHamiltonian, {TXASyz,TXASxz,TXASxy,TXASz2,TXASx2y2}, psiList, {{"Emin",-5}, {"Emax",10}, {"NE",3000}, {"Gamma",0.1}})
-- 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
XASIsoSpectra = Spectra.Sum(XASSpectra,{1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,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
XASIsoSpectra.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.
XASSpectra.Print({{"file","XASSpec.dat"}})
XASIsoSpectra.Print({{"file","XASIsoSpec.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 'XASSpec.ps'
set size 1.0, 1.0
set terminal postscript landscape enhanced color  "Times" 8

set xrange [8332.5:8335]

energyshift=8331

set multiplot layout 5, 3

plot "XASSpec.dat"  u ($1+energyshift):(-$3 ) title '  yz-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$5 ) title '  yz-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$7 ) title '  yz-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$9 ) title '  xz-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$11) title '  xz-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$13) title '  xz-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$15) title '  xy-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$17) title '  xy-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$19) title '  xy-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$21) title 'x2y2-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$23) title 'x2y2-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$25) title 'x2y2-polarized Sz= 1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$27) title '  z2-polarized Sz=-1' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$29) title '  z2-polarized Sz= 0' with lines ls  1
plot "XASSpec.dat"  u ($1+energyshift):(-$31) title '  z2-polarized Sz= 1' with lines ls  1

unset multiplot

set xrange [8330:8336]

intensityscale=64

plot "XASSpec.dat" using ($1+energyshift):((-$3-$9-$15-$21-$27) * intensityscale) title 'isotropic theory' with lines ls 1


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

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

]]

-- next we save this string 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 I transform the format
os.execute(" ps2pdf XASSpec.ps ; ps2pdf XASIsoSpec.ps")