-- For NiO we know that 10Dq=1.1 eV, however if you have a new compound you might roughly
-- know these values, but never exactly. It is then nice to see how the eigen-state
-- energies change as a function of 10Dq.

-- plots were one shows the eigen-state energy as a function of some parameter (crystal-
-- field strength) are called energy level diagrams or Sugano - Tanabe diagrams 

-- here we create the diagram for Ni - d8

-- in order to minimize the output of the calculation we set the verbosity to 0
Verbosity(0)

-- we have a single d-shell in the bases and thus need 10 Fermions, grouped in 5 orbitals
-- with spin down and 5 with spin up.
NF=10
NB=0
IndexDn_3d={0,2,4,6,8}
IndexUp_3d={1,3,5,7,9}

-- again we define several operators

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)

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

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

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)

-- for d^8 there are 10*9/2 = 45 different states. For a diagram we want to calculate all
Npsi=45

-- as before we can set the parameters, in this case all but 10Dq
U       =  0.000
F2dd    = 11.142 
F4dd    =  6.874
F0dd    = U+(F2dd+F4dd)*2/63
zeta_3d = 0.081
Bz      = 0.000001

-- and we define the Hamiltonian as a sum of prefactors times operators. Note that
-- we do not include the crystal field here.
Hamiltonian0 =  F0dd*OppF0 + F2dd*OppF2 + F4dd*OppF4 + zeta_3d*Oppldots + Bz*(2*OppSz + OppLz)

-- We first calculate the eigen states of Hamiltonian0. For a filling of 8 electrons
-- the wavefunctions are stored as a table in the variable psiList
StartRestrictions = {NF, NB,  {"1111111111",8,8}}
psiList = Eigensystem(Hamiltonian0, StartRestrictions, Npsi)

-- We will write the eigen energies to the file EnergyLevelDiagram, which we open for writing
-- The variable file now points to the file "EnergyLevelDiagram"
file = assert(io.open("EnergyLevelDiagram", "w"))

-- no we need to loop over the value of 10Dq and for each value diagonalize the Hamiltonian
-- and save the eigenvalues to the file EnergyLevelDiagram
-- we do 301 calculations changing 10Dq from 0 to 3 in steps of 10 meV.
-- i is an index that counts from 0 to 300
-- tenDq is a variable equal to 0.01*i and the actual crystal-field strength.
-- During the loop we write the value of 10Dq to the screen so we know where the calculation is
-- Into the file we write first the value of 10Dq and then the 45 eigenvalues
-- but in order to calculate the eigenvalues we first need to call the function Eigensystem
-- acting on Hamiltonian0 + tenDq * OpptenDq
-- Eigensystem is a flexible function that can be used in different ways and on different
-- objects. If you have a set of states close to the eigenstates you are interested in you can
-- use these states as starting point.
-- Calling Eigensystem with an operator as the first argument and a table of functions as the second
-- will change the functions to become the lowest eigenstates of the operator.
io.write("10Dq = ")
for i=0, 300 do
  tenDq = 0.01*i
  io.write(string.format("%5.3f ",tenDq))
  io.flush()
  file:write(string.format("%14.7E ",tenDq))
  Hamiltonian=Hamiltonian0 + tenDq * OpptenDq
  Eigensystem(Hamiltonian, psiList)
  for key,value in pairs(psiList) do
    energy = value * Hamiltonian * value
    file:write(string.format("%14.7E ",energy))
  end
  file:write("\n")
end
io.write("\n")
file:close()

-- Once we have created a file containing the eigenenergies we can make a plot
-- A gnuplot script that would make this plot can be found below
gnuplotInput = [[
set autoscale 
set xtic auto  
set ytic auto  
set style line  1 lt 1 lw 1 lc rgb "#000000"

set xlabel "10Dq (eV)" font "Times,12"
set ylabel "Energy (eV)" font "Times,12"

set out 'EnergyLevelDiagram.ps'
set size 1.0, 0.625
set terminal postscript portrait enhanced color  "Times" 12

plot for [i=2:46] "EnergyLevelDiagram" using 1:i notitle with lines ls  1
]]

-- write the gnuplot script to a file
file = io.open("EnergyLevelDiagram.gnuplot", "w")
file:write(gnuplotInput)
file:close()

-- call gnuplot to execute the script
os.execute("gnuplot EnergyLevelDiagram.gnuplot")
-- change the postscript file to pdf or eps
os.execute("ps2pdf EnergyLevelDiagram.ps ; ps2eps EnergyLevelDiagram.ps ; mv EnergyLevelDiagram.eps temp.eps ; eps2eps temp.eps EnergyLevelDiagram.eps ; rm temp.eps")