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 molecule\n")
TB = TightBindingDefFromDresdenFPLO(FPLOOut)

-- For the CrF4 cluster we can rotate the ligand orbitals
print("--Rotate the Ligand orbitals--\n")
RTB, T = BlockBandDiagonalize(TB,{{"Cr",{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}})
--give the indices in the system a name and calculate the number of fermionic states:
Index, NFermi = CreateAtomicIndicesDict({"Cr_3d","Ligand_d","Ligand_p"},{{"Ligand",{"Ligand_d","Ligand_p"}}})


-- 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
rhoNoSpind = Matrix.Sub(rhoNoSpin,5)
-- add spin
rhod = Matrix.AddSpin(rhoNoSpind)

-- 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["Cr_3d_Up"],Index["Cr_3d_Dn"],{1,0,0}),YtoZMatrix("d"))
OppF2 = Rotate(NewOperator("U",10,Index["Cr_3d_Up"],Index["Cr_3d_Dn"],{0,1,0}),YtoZMatrix("d"))
OppF4 = Rotate(NewOperator("U",10,Index["Cr_3d_Up"],Index["Cr_3d_Dn"],{0,0,1}),YtoZMatrix("d"))

-- 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,rhod,{{"AddDFTSelfInteraction",true}}))
OppF2MFDFT = Chop(MeanFieldOperator(OppF2,rhod,{{"AddDFTSelfInteraction",true}}))
OppF4MFDFT = Chop(MeanFieldOperator(OppF4,rhod,{{"AddDFTSelfInteraction",true}}))

-- change the dimension of the Coulomb operator the work on a space that includes both the d and the ligand orbitals
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 = 2
Delta   =  2
Udd     =  4.5
--The slater integrals are calculated using the radial functions generated by FPLO, which can be read within Quanty:
correlatedRadialFunctions = ReadFPLOBasisFunctions({"3d"},"../DFT/+fval.001.1")
slaterIntegrals = GetSlaterIntegrals({"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
-- Magnetic field
Bz = 0.1 * EnergyUnits.Tesla.value
-- Spin orbit
zeta_3d =  0.041

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

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

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

-- define some aditional operators
print("--Define some additional operators--\n")
YtoZdouble = YtoZMatrix({"3d","Ligand_d","Ligand_p"})

OppSz    = Rotate(NewOperator("Sz",   NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)
OppLz    = Rotate(NewOperator("Lz",   NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)
Oppldots = Rotate(NewOperator("ldots",NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)

OppSsqr  = Rotate(NewOperator("Ssqr", NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)
OppLsqr  = Rotate(NewOperator("Lsqr", NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)
OppJsqr  = Rotate(NewOperator("Jsqr", NFermi,Index["Cr_3d_Up"],Index["Cr_3d_Dn"]),YtoZdouble)

OppN_Ni      = NewOperator("Number",NFermi,Index["Cr_3d"],                Index["Cr_3d"],                {1,1,1,1,1,1,1,1,1,1})
OppN_Ligand  = NewOperator("Number",NFermi,Index["Ligand"],               Index["Ligand"],               {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1})
OppN_Ni_yz   = NewOperator("Number",NFermi,IndexFPLO["0_Cr_3d_d_{yz}"],       IndexFPLO["0_Cr_3d_d_{yz}"],       {1,1})
OppN_Ni_xz   = NewOperator("Number",NFermi,IndexFPLO["0_Cr_3d_d_{xz}"],       IndexFPLO["0_Cr_3d_d_{xz}"],       {1,1})
OppN_Ni_xy   = NewOperator("Number",NFermi,IndexFPLO["0_Cr_3d_d_{xy}"],       IndexFPLO["0_Cr_3d_d_{xy}"],       {1,1})
OppN_Ni_x2y2 = NewOperator("Number",NFermi,IndexFPLO["0_Cr_3d_d_{x^2-y^2}"],  IndexFPLO["0_Cr_3d_d_{x^2-y^2}"],  {1,1})
OppN_Ni_z2   = NewOperator("Number",NFermi,IndexFPLO["0_Cr_3d_d_{3z^2-r^2}"], IndexFPLO["0_Cr_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

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

-- correct for screening: set the onsite energy of the Cr_3d and Ligand shell to ed and eL
print("Correct for screening: set the onsite energy of the Cr_3d and Ligand shell to ed and eL\n")
ed      = (16*Delta-nd*(31+nd)*Udd/2)/(16+nd)
eL      = nd*((1+nd)*Udd/2-Delta)/(16+nd)
print("ed = "..ed.." eV")
print("eL = "..eL.." eV")
OperatorSetTrace(Hamiltonian,ed,Index["Cr_3d"])
OperatorSetTrace(Hamiltonian,eL,Index["Ligand"])


print("--Compute eigenstates--\n")
-- we now can create the lowest Npsi eigenstates:
Npsi=30
-- in order to make sure we have a filling of 8 electrons we need to define some restrictions
StartRestrictions = {NFermi, 0, {DeterminantString(NFermi,Index["Cr_3d"]),nd,nd}, {DeterminantString(NFermi,Index["Ligand"]),16,16}}
print("Start Restrictions:")
print(StartRestrictions)
Restrictions = {NFermi, 0, {DeterminantString(NFermi,Index["Cr_3d"]),nd,nd+2}, {DeterminantString(NFermi,Index["Ligand"]),16-2,16}}
print("Restrictions:")
print(Restrictions)

psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi,{{"Restrictions", Restrictions}})

-- 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.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, OppN_Ni_eg, OppN_Ni_t2g, OppN_Ni, OppN_Ligand}

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