
dofile("../Expand_Ligand_Basis.Quanty")

Verbosity(0)
-- This tutorial calculates the ground-state of NiO within the Anderson impurity approximation.
-- We start with a ligand field theory model and then replace the single bath site by a band
-- Have a look at the ligand field theory calculations first.

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

ed      = (10*Delta-nd*(19+nd)*U/2)/(10+nd)
eL      = nd*((1+nd)*U/2-Delta)/(10+nd)
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 = 8
-- 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 = 0

-- derived parameters
eL           = eL + eLshift
epsilon0Bath = eL * Matrix.Identity(5)
tMat         = (W/4) * Matrix.Identity(5)

EshiftTotal = 10*NShells*eL

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

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

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

-- Coulomb interaction

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

-- crystal and ligand field

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)

-- Neg and Nt2g operators

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
            + ed * OppN_3d + eL * 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

-- in orde to see the effect of the bath we can look at the density matrix
print("The digaonal part of the one particle density matrix")
rho = DensityMatrix(psiList[1])
for i=1, #rho do
  print(string.format("%3d ",i),rho[i][i])
end