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