Verbosity(0)

-- read the output of FPLO
print("--Read FPLO output--\n")
FPLOOut = FileReadDresdenFPLO("../DFT/out.wan")

-- from the DFT output we can create a tight binding Hamiltonian
print("--Create the Tight Binding Hamiltonian--\n")
print("Create the tight binding Hamiltonian for the crystal\n")
TB = TightBindingDefFromDresdenFPLO(FPLOOut)

print("Define a cluster")
NewCluster = FindAllAtomsInsideSphere(TB.Atoms,TB.Cell,{0,0,0},4)
print("Quanty detected the following cluster:")
print(NewCluster)
print("")

-- from the TB output of a crystal we can create a TB object of a Cluster
print("Create the tight binding Hamiltonian for the cluster\n")

HDFTLarge, ClusterTB, IndexName = CreateClusterHamiltonian(TB, {"open", NewCluster},  {{"AddSpin",false}, {"ReturnTBIndicesDict",true}})

-- For the NiO6 cluster we can rotate the ligand orbitals
print("--Rotate the Ligand orbitals--\n")
RTB, T = BlockBandDiagonalize(ClusterTB,{{"Ni",{0,0,0}}})

-- from the rotated TB output we can create a Cluster Hamiltonian (Operator instead of TB object)
print("Create Cluster Hamiltonian after rotation\n")
HDFT, ClusterTB, IndexFPLO = CreateClusterHamiltonian(RTB, {"open", RTB.Atoms},  {{"AddSpin",true}, {"ReturnTBIndicesDict",true}})

--Now we want to add the Ni 2_p shell to the system
--New index object (first argument is a list of all orbitals, the second argument groups the two Ni orbitals)
Index, NFermi = CreateAtomicIndicesDict({"Ni_3d","Ligand_d","Ni_2p"}, {{"Ni",{"Ni_2p","Ni_3d"}}})
HDFT.NF = NFermi
print("The indices of the system are now fixed:")
print(Index)
print("")

-- Next we need to correct for the double counting (i.e. subtract the DFT Mean-field approximation of the Coulomb interaction)
print("--Subtract the DFT Mean-field Coulomb interaction between the d electrons from the tight-binding Hamiltonian--\n")
-- calculate DFT density Matrix in order to create double counting potential
print("Calculate the DFT 1-particle density matrix\n")
rhoNoSpin = Chop(CalculateRho(TB))
-- reduce the matrix to keep only the d shell
rhoNoSpin_d = Matrix.ToTable(Matrix.Sub(rhoNoSpin,5))
-- add spin
rho_d = Matrix.AddSpin(rhoNoSpin_d)

-- create the Coulomb operator on a basis of tesseral Harmonics
print("Create the Coulomb operator for the d-shell\n")
OppF0 = Rotate(NewOperator("U",10,Index["Ni_3d_Up"],Index["Ni_3d_Dn"],{1,0,0}),YtoZMatrix({"d"}))
OppF2 = Rotate(NewOperator("U",10,Index["Ni_3d_Up"],Index["Ni_3d_Dn"],{0,1,0}),YtoZMatrix({"d"}))
OppF4 = Rotate(NewOperator("U",10,Index["Ni_3d_Up"],Index["Ni_3d_Dn"],{0,0,1}),YtoZMatrix({"d"}))

print("Create the Coulomb operator between the p- and d-shell\n")
OppUpdF0 = Rotate(NewOperator("U", NFermi, Index["Ni_2p_Up"],Index["Ni_2p_Dn"], Index["Ni_3d_Up"],Index["Ni_3d_Dn"], {1,0}, {0,0}),YtoZMatrix({"Ni_3d","Ligand_d","Ni_2p"}))
OppUpdF2 = Rotate(NewOperator("U", NFermi, Index["Ni_2p_Up"],Index["Ni_2p_Dn"], Index["Ni_3d_Up"],Index["Ni_3d_Dn"], {0,1}, {0,0}),YtoZMatrix({"Ni_3d","Ligand_d","Ni_2p"}))
OppUpdG1 = Rotate(NewOperator("U", NFermi, Index["Ni_2p_Up"],Index["Ni_2p_Dn"], Index["Ni_3d_Up"],Index["Ni_3d_Dn"], {0,0}, {1,0}),YtoZMatrix({"Ni_3d","Ligand_d","Ni_2p"}))
OppUpdG3 = Rotate(NewOperator("U", NFermi, Index["Ni_2p_Up"],Index["Ni_2p_Dn"], Index["Ni_3d_Up"],Index["Ni_3d_Dn"], {0,0}, {0,1}),YtoZMatrix({"Ni_3d","Ligand_d","Ni_2p"}))

-- create the Hartree-Fock mean-field approximation that needs to be substracted from the DFT results
print("Create the Coulomb operator for the d-shell using the mean-field approximation (needs the DFT density)\n")
OppF0MFDFT = Chop(MeanFieldOperator(OppF0,rho_d,{{"AddDFTSelfInteraction",true}}))
OppF2MFDFT = Chop(MeanFieldOperator(OppF2,rho_d,{{"AddDFTSelfInteraction",true}}))
OppF4MFDFT = Chop(MeanFieldOperator(OppF4,rho_d,{{"AddDFTSelfInteraction",true}}))

-- change the dimension of the Coulomb operator the work on a space that includes both the d and the ligand orbitals
--The operators acting on the d- and p-shell furthermore need to be rotated
OppF0.NF = NFermi
OppF2.NF = NFermi
OppF4.NF = NFermi
OppF0MFDFT.NF = NFermi
OppF2MFDFT.NF = NFermi
OppF4MFDFT.NF = NFermi

-- set parameters
print("--Set the parameters--\n")
nd = 8
Delta   =  4.7
Udd       =  7.30
Upd     =  8.5
--The slater integrals are calculated using the radial functions generated by FPLO, which can be read within Quanty:
correlatedRadialFunctionsp = ReadFPLOBasisFunctions({"2p"},"../DFT/+fcor.001.1")
correlatedRadialFunctionsd = ReadFPLOBasisFunctions({"3d"},"../DFT/+fval.001.1")
--The radial functions are saved in two different files, but we need them in the same object:
correlatedRadialFunctions={correlatedRadialFunctionsp[1],correlatedRadialFunctionsd[1]}
--Now the slater integrals can be calculated
slaterIntegrals = GetSlaterIntegrals({"2p","3d"},correlatedRadialFunctions)

--The resulting object can then be evaluated at the position of interest:
F2dd = slaterIntegrals["3d 3d 3d 3d"][2] * EnergyUnits.Ha.value
F4dd = slaterIntegrals["3d 3d 3d 3d"][4] * EnergyUnits.Ha.value
F0dd    = Udd+(F2dd+F4dd)*2/63
F2pd    =  slaterIntegrals["3d 2p 3d 2p"][2] * EnergyUnits.Ha.value
G1pd    =  slaterIntegrals["2p 3d 3d 2p"][1] * EnergyUnits.Ha.value
G3pd    =  slaterIntegrals["2p 3d 3d 2p"][3] * EnergyUnits.Ha.value
F0pd    = Upd + (1/15)*G1pd + (3/70)*G3pd

-- Magnetic field
Bz = 0.1 * EnergyUnits.Tesla.value
Hz      =  0.120
-- Spin orbit
zeta_3d =  0.081
zeta_2p = 11.51

print("Parameters:")
print("nd = "..nd)
print("Delta = "..Delta.." eV")
print("Udd = "..Udd.." eV")
print("Upd = "..Upd.." eV")
print("Bz = "..Bz.." eV")
print("Hz = "..Hz.." eV")
print("zeta_3d = "..zeta_3d.." eV")
print("zeta_2p = "..zeta_2p.." eV")

print("Slater integrals computed from DFT radial wavefunctions:")
print("F2dd = "..F2dd.." eV")
print("F4dd = "..F4dd.." eV")
print("F2pd = "..F2pd.." eV")
print("G1pd = "..G1pd.." eV")
print("G3pd = "..G3pd.." eV")

print("Monopole like Slater integrals (from U):")
print("F0dd = "..F0dd.." eV")
print("F0pd = "..F0pd.." eV")
print("")



-- define some aditional operators
print("--Define some additional operators--\n")
YtoZtriple = YtoZMatrix({"Ni_3d","Ligand_d","Ni_2p"})

OppSz    = Rotate(NewOperator("Sz",   NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)
OppLz    = Rotate(NewOperator("Lz",   NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)
Oppldots_3d = Rotate(NewOperator("ldots",NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)
Oppldots_2p = Rotate(NewOperator("ldots",NFermi,Index["Ni_2p_Up"],Index["Ni_2p_Dn"]),YtoZtriple)

OppSsqr  = Rotate(NewOperator("Ssqr", NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)
OppLsqr  = Rotate(NewOperator("Lsqr", NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)
OppJsqr  = Rotate(NewOperator("Jsqr", NFermi,Index["Ni_3d_Up"],Index["Ni_3d_Dn"]),YtoZtriple)

OppN_Ni      = NewOperator("Number",NFermi,Index["Ni"],                       Index["Ni"],                       {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1})
OppN_Ligand  = NewOperator("Number",NFermi,Index["Ligand_d"],                 Index["Ligand_d"],                 {1,1,1,1,1,1,1,1,1,1})
OppN_Ni_yz   = NewOperator("Number",NFermi,IndexFPLO["0_Ni_3d_d_{yz}"],       IndexFPLO["0_Ni_3d_d_{yz}"],       {1,1})
OppN_Ni_xz   = NewOperator("Number",NFermi,IndexFPLO["0_Ni_3d_d_{xz}"],       IndexFPLO["0_Ni_3d_d_{xz}"],       {1,1})
OppN_Ni_xy   = NewOperator("Number",NFermi,IndexFPLO["0_Ni_3d_d_{xy}"],       IndexFPLO["0_Ni_3d_d_{xy}"],       {1,1})
OppN_Ni_x2y2 = NewOperator("Number",NFermi,IndexFPLO["0_Ni_3d_d_{x^2-y^2}"],  IndexFPLO["0_Ni_3d_d_{x^2-y^2}"],  {1,1})
OppN_Ni_z2   = NewOperator("Number",NFermi,IndexFPLO["0_Ni_3d_d_{3z^2-r^2}"], IndexFPLO["0_Ni_3d_d_{3z^2-r^2}"], {1,1})

OppN_Ni_eg  = OppN_Ni_x2y2 + OppN_Ni_z2
OppN_Ni_t2g = OppN_Ni_yz + OppN_Ni_xz + OppN_Ni_xy

-- dipole transition
 
t=math.sqrt(1/2)
 
Akm = {{1,-1,t},{1, 1,-t}}
TXASx = NewOperator("CF", NFermi, Index["Ni_3d_Up"], Index["Ni_3d_Dn"], Index["Ni_2p_Up"], Index["Ni_2p_Dn"], Akm)
Akm = {{1,-1,t*I},{1, 1,t*I}}
TXASy = NewOperator("CF", NFermi, Index["Ni_3d_Up"], Index["Ni_3d_Dn"], Index["Ni_2p_Up"], Index["Ni_2p_Dn"], Akm)
Akm = {{1,0,1}}
TXASz = NewOperator("CF", NFermi, Index["Ni_3d_Up"], Index["Ni_3d_Dn"], Index["Ni_2p_Up"], Index["Ni_2p_Dn"], Akm)
 
TXASr = t*(TXASx - I * TXASy)
TXASl =-t*(TXASx + I * TXASy)

-- define the Hamiltonian
print("--Define Hamiltonian--\n")
Hamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
            + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
            + Bz * (2*OppSz + OppLz) + zeta_3d * Oppldots_3d + Hz * OppSz
            
-- correct for screening: set the onsite energy of the Ni_3d and Ligand shell to ed and eL
print("Correct for screening: set the onsite energy of the Ni_3d and Ligand shell to ed and eL\n")
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)
print("ed = "..ed.." eV")
print("eL = "..eL.." eV")
print("epfinal = "..epfinal.." eV")
print("edfinal = "..edfinal.." eV")
print("eLfinal = "..eLfinal.." eV")
OperatorSetTrace(Hamiltonian,0,Index["Ni_2p"])
OperatorSetTrace(Hamiltonian,ed,Index["Ni_3d"])
OperatorSetTrace(Hamiltonian,eL,Index["Ligand_d"])
            
print("--Define XAS-Hamiltonian--\n")
XASHamiltonian = HDFT - F0dd * OppF0MFDFT - F2dd * OppF2MFDFT - F4dd * OppF4MFDFT
            + F0dd * OppF0 + F2dd * OppF2 + F4dd * OppF4
            + Bz * (2*OppSz + OppLz) + zeta_3d * Oppldots_3d + Hz * OppSz
            + zeta_2p * Oppldots_2p
            + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3
            

print("Correct for screening: set the onsite energy of the Ni_2p, Ni_3d and Ligand shell to epfinal, edfinal and eLfinal\n")
OperatorSetTrace(XASHamiltonian,epfinal,Index["Ni_2p"])
OperatorSetTrace(XASHamiltonian,edfinal,Index["Ni_3d"])
OperatorSetTrace(XASHamiltonian,eLfinal,Index["Ligand_d"])


print("--Compute eigenstates--\n")
-- we now can create the lowest Npsi eigenstates:
Npsi=3
-- in order to make sure we have a filling of 8 electrons we need to define some restrictions
StartRestrictions = {NFermi, 0, {DeterminantString(NFermi,Index["Ni_3d"]),nd,nd}, {DeterminantString(NFermi,Index["Ligand_d"]),10,10}, {DeterminantString(NFermi,Index["Ni_2p"]),6,6}}
print("Start Restrictions:")
print(StartRestrictions)

psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)

-- print of some expectation values
Hamiltonian.Name = "<E>"
OppSsqr.Name = "<S^2>"
OppLsqr.Name = "<L^2>"
OppJsqr.Name = "<J^2>"
OppSz.Name = "<S_z^3d>"
OppLz.Name = "<L_z^3d>"
Oppldots_3d.Name = "<l.s>"
OppN_Ni_eg.Name = "<Neg_Ni>"
OppN_Ni_t2g.Name = "<Nt2g_Ni>"
OppN_Ni.Name = "<N_Ni>"
OppN_Ligand.Name = "<N_L>"
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, Oppldots_3d, OppN_Ni_eg, OppN_Ni_t2g, OppN_Ni, OppN_Ligand}

PrintExpectationValues(psiList,oppList,{{"ColWidth",10}})


print("--Create the Spectra--\n")
-- 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}})
 
-- 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"}})
 
 
print("--Plot the Spectra--\n")
-- 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.6
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 ; ps2eps XASIsoSpec.ps ;  mv XASIsoSpec.eps temp.eps ; eps2eps temp.eps XASIsoSpec.eps ; rm temp.eps")
