Differences

This shows you the differences between two versions of the page.

Link to this comparison view

documentation:tutorials:nio_crystal_field:energy_level_diagram [2016/10/08 19:11] – created Maurits W. Haverkortdocumentation:tutorials:nio_crystal_field:energy_level_diagram [2016/10/10 09:41] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +{{indexmenu_n>3}}
 +====== Energy level diagram ======
  
 +###
 +In order to do temperature averaging it is important to understand the number of excited states that are important. One can learn a lot by looking at the energy level diagram. Here we plot one for Ni$^{2+}$.
 +###
 +
 +###
 +The input file is:
 +<code Quanty Energy_level_diagram.Quanty>
 +-- 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")
 +</code>
 +###
 +
 +###
 +As in the [[documentation:tutorials:small_programs_a_quick_start:energy_level_diagram|example in the beginning of the tutorial]] Quanty returns a nice plot. Note that one can add labeling. For this have a look at the [[documentation:tutorials:small_programs_a_quick_start:energy_level_diagram|example in the beginning of the tutorial]].
 +
 +{{:documentation:tutorials:nio_crystal_field:energyleveldiagram.png?nolink |}}
 +###
 +
 +###
 +The text output is trivial and just shows counters that tell you in which loop you are during the calculation:
 +<file Quanty_Output Energy_level_diagram.out>
 +10Dq = 0.000 0.010 0.020 0.030 0.040 0.050 0.060 0.070 0.080 0.090 0.100 0.110 0.120 0.130 0.140 0.150 0.160 0.170 0.180 0.190 0.200 0.210 0.220 0.230 0.240 0.250 0.260 0.270 0.280 0.290 0.300 0.310 0.320 0.330 0.340 0.350 0.360 0.370 0.380 0.390 0.400 0.410 0.420 0.430 0.440 0.450 0.460 0.470 0.480 0.490 0.500 0.510 0.520 0.530 0.540 0.550 0.560 0.570 0.580 0.590 0.600 0.610 0.620 0.630 0.640 0.650 0.660 0.670 0.680 0.690 0.700 0.710 0.720 0.730 0.740 0.750 0.760 0.770 0.780 0.790 0.800 0.810 0.820 0.830 0.840 0.850 0.860 0.870 0.880 0.890 0.900 0.910 0.920 0.930 0.940 0.950 0.960 0.970 0.980 0.990 1.000 1.010 1.020 1.030 1.040 1.050 1.060 1.070 1.080 1.090 1.100 1.110 1.120 1.130 1.140 1.150 1.160 1.170 1.180 1.190 1.200 1.210 1.220 1.230 1.240 1.250 1.260 1.270 1.280 1.290 1.300 1.310 1.320 1.330 1.340 1.350 1.360 1.370 1.380 1.390 1.400 1.410 1.420 1.430 1.440 1.450 1.460 1.470 1.480 1.490 1.500 1.510 1.520 1.530 1.540 1.550 1.560 1.570 1.580 1.590 1.600 1.610 1.620 1.630 1.640 1.650 1.660 1.670 1.680 1.690 1.700 1.710 1.720 1.730 1.740 1.750 1.760 1.770 1.780 1.790 1.800 1.810 1.820 1.830 1.840 1.850 1.860 1.870 1.880 1.890 1.900 1.910 1.920 1.930 1.940 1.950 1.960 1.970 1.980 1.990 2.000 2.010 2.020 2.030 2.040 2.050 2.060 2.070 2.080 2.090 2.100 2.110 2.120 2.130 2.140 2.150 2.160 2.170 2.180 2.190 2.200 2.210 2.220 2.230 2.240 2.250 2.260 2.270 2.280 2.290 2.300 2.310 2.320 2.330 2.340 2.350 2.360 2.370 2.380 2.390 2.400 2.410 2.420 2.430 2.440 2.450 2.460 2.470 2.480 2.490 2.500 2.510 2.520 2.530 2.540 2.550 2.560 2.570 2.580 2.590 2.600 2.610 2.620 2.630 2.640 2.650 2.660 2.670 2.680 2.690 2.700 2.710 2.720 2.730 2.740 2.750 2.760 2.770 2.780 2.790 2.800 2.810 2.820 2.830 2.840 2.850 2.860 2.870 2.880 2.890 2.900 2.910 2.920 2.930 2.940 2.950 2.960 2.970 2.980 2.990 3.000 
 +</file>
 +###
 +
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » NiO crystal field » Energy level diagram

Print/export