Differences

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

Link to this comparison view

documentation:tutorials:nio_ligand_field:xas_l23_partial_excitations [2016/10/09 15:49] – created Maurits W. Haverkortdocumentation:tutorials:nio_ligand_field:xas_l23_partial_excitations [2016/10/10 09:41] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +{{indexmenu_n>14}}
 +====== XAS $L_{2,3}$ partial excitations ======
  
 +###
 +In order to understand where a particular peak originates from one can look at the character of the excited state. Quanty works with Green's functions and we can not just look at the excited eigen-states as these are not really calculated. We can however modify the transition operator in order to understand where the spectra come from. The example below looks at excitations into the $t_{2g}$ and $e_g$ irreducible representations of the $d$ orbitals separately and looks at excitations from the $j=1/2$ and $j=3/2$ sub-shells of the $2p$ orbitals separately. It furthermore looks at excitations into the $d^9$ and $d^{10}$ configurations.
 +###
 +
 +###
 +The input file is:
 +<code Quanty XAS_partial.Quanty>
 +-- In order to understand the different peaks in an absorption experiment one can look
 +-- at the character of the states one excites into. 
 +
 +-- here we calculate the 2p to 3d x-ray absorption of NiO within the Ligand-field theory
 +-- approximation. The script is basically a copy of tutorial 41. The difference is in how
 +-- we define the transition operator.
 +
 +-- we calculate the spectra separated to t2g and eg excitations, to 2p-j=1/2 or 2p-j=3/2
 +-- excitations and to excitations into d9 L^10 or d10L9 states.
 +
 +-- from the previous example we know that within NiO there are 3 states close to each other
 +-- and then there is an energy gap of about 1 eV. We thus only need to consider the 3
 +-- lowest states (Npsi=3 later on)
 +
 +NF=26
 +NB=0
 +IndexDn_2p={ 0, 2, 4}
 +IndexUp_2p={ 1, 3, 5}
 +IndexDn_3d={ 6, 8,10,12,14}
 +IndexUp_3d={ 7, 9,11,13,15}
 +IndexDn_Ld={16,18,20,22,24}
 +IndexUp_Ld={17,19,21,23,25}
 +
 +-- angular momentum operators on the d-shell
 +
 +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)
 +
 +-- Angular momentum operators on the Ligand 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)
 +
 +-- total angular momentum
 +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
 +
 +-- define the coulomb operator
 +-- we here define the part depending on F0 seperately from the part depending on F2
 +-- when summing we can put in the numerical values of the slater integrals
 +
 +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})
 +
 +-- define onsite energies - crystal field
 +-- Akm = {{k1,m1,Akm1},{k2,m2,Akm2}, ... }
 +
 +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)
 +
 +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)
 +
 +OppNUp_2p = NewOperator("Number", NF, IndexUp_2p, IndexUp_2p, {1,1,1})
 +OppNDn_2p = NewOperator("Number", NF, IndexDn_2p, IndexDn_2p, {1,1,1})
 +OppN_2p = OppNUp_2p + OppNDn_2p
 +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
 +
 +-- define L-d interaction
 +
 +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)
 +
 +-- core valence interaction
 +
 +Oppcldots= NewOperator("ldots", NF, IndexUp_2p, IndexDn_2p)
 +OppUpdF0 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {1,0}, {0,0})
 +OppUpdF2 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,1}, {0,0})
 +OppUpdG1 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {1,0})
 +OppUpdG3 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {0,1})
 +
 +-- dipole transition
 +
 +t=math.sqrt(1/2)
 +
 +Akm = {{1,-1,t},{1, 1,-t}}
 +TXASx = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 +Akm = {{1,-1,t*I},{1, 1,t*I}}
 +TXASy = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 +Akm = {{1,0,1}}
 +TXASz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 +
 +-- we now define rotation matrices that rotate the basis of (l lz s sz) to a basis of
 +-- (k tau s sz) for the d-shell or (j jz) for the p-shell. The function rotate calculates
 +-- R^* Opp R^T and the rotation matrix is thus equal to the eigen-functions of the
 +-- crystal-field operator and the spin-orbit coupling operator respectively.
 +
 +-- the total matrix is a 16 by 16 matrix
 +
 +t=sqrt(1/2)
 +u=I*sqrt(1/2)
 +d=sqrt(1/3)
 +e=sqrt(2/3)
 +
 +-- rotating the 2p states to a (j jz) basis
 +
 +rotmatjjzp = {{ 0,-e, d, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0,-d, e, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, d, e, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, e, d, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +              { 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +              
 +-- rotating the 3d states to a (k tau s sz) basis
 +              
 +rotmatKd   = {{ 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +              { 0, 0, 0, 0, 0, 0,     t, 0, 0, 0, 0, 0, 0, 0, t, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, t, 0, 0, 0, 0, 0, 0, 0, t,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, u, 0, 0, 0, u, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, u, 0, 0, 0, u, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, t, 0, 0, 0,-t, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, t, 0, 0, 0,-t, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     u, 0, 0, 0, 0, 0, 0, 0,-u, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, u, 0, 0, 0, 0, 0, 0, 0,-u,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +              { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +              
 +-- We can modify the rotation of the 2p orbitals such that we only
 +-- keep the j=1/2 sub-shell
 +
 +rotmatjjzp12 = {{ 0,-e, d, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0,-d, e, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +                
 +-- or modify it such that we only keep the j=3/2 sub-shell
 +                
 +rotmatjjzp32 = {{ 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, d, e, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, e, d, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +                
 +-- We can modify the rotation of the 3d orbitals such that we only keep
 +-- the orbitals belonging to the eg irreducible representation.
 +                
 +rotmatKdeg   = {{ 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     t, 0, 0, 0, 0, 0, 0, 0, t, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, t, 0, 0, 0, 0, 0, 0, 0, t,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +                
 +-- or such that we only keep the orbitals belonging to the t2g irreducible representation
 +                
 +rotmatKdt2g  = {{ 1, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 1, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 1, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 1, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, u, 0, 0, 0, u, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, u, 0, 0, 0, u, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, t, 0, 0, 0,-t, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, t, 0, 0, 0,-t, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     u, 0, 0, 0, 0, 0, 0, 0,-u, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, u, 0, 0, 0, 0, 0, 0, 0,-u,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +              
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 1, 0, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
 +                { 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}
 +                
 +-- The transition operator restricted to a subset of orbitals is then given as:
 +-- Tnew = R2^* . R1^* . Told . R1^T . R2^T
 +-- with R1 the restricted rotation matrix
 +-- and R2 the full rotation matrix conjugate transposed (R^+)
 +
 +TXASxj12 = Chop(Rotate(Chop(Rotate(TXASx,rotmatjjzp12)),ConjugateTranspose(rotmatjjzp)))
 +TXASyj12 = Chop(Rotate(Chop(Rotate(TXASy,rotmatjjzp12)),ConjugateTranspose(rotmatjjzp)))
 +TXASzj12 = Chop(Rotate(Chop(Rotate(TXASz,rotmatjjzp12)),ConjugateTranspose(rotmatjjzp)))
 +TXASxj32 = Chop(Rotate(Chop(Rotate(TXASx,rotmatjjzp32)),ConjugateTranspose(rotmatjjzp)))
 +TXASyj32 = Chop(Rotate(Chop(Rotate(TXASy,rotmatjjzp32)),ConjugateTranspose(rotmatjjzp)))
 +TXASzj32 = Chop(Rotate(Chop(Rotate(TXASz,rotmatjjzp32)),ConjugateTranspose(rotmatjjzp)))
 +TXASxeg  = Chop(Rotate(Chop(Rotate(TXASx,rotmatKdeg  )),ConjugateTranspose(rotmatKd  )))
 +TXASyeg  = Chop(Rotate(Chop(Rotate(TXASy,rotmatKdeg  )),ConjugateTranspose(rotmatKd  )))
 +TXASzeg  = Chop(Rotate(Chop(Rotate(TXASz,rotmatKdeg  )),ConjugateTranspose(rotmatKd  )))
 +TXASxt2g = Chop(Rotate(Chop(Rotate(TXASx,rotmatKdt2g )),ConjugateTranspose(rotmatKd  )))
 +TXASyt2g = Chop(Rotate(Chop(Rotate(TXASy,rotmatKdt2g )),ConjugateTranspose(rotmatKd  )))
 +TXASzt2g = Chop(Rotate(Chop(Rotate(TXASz,rotmatKdt2g )),ConjugateTranspose(rotmatKd  )))
 +
 +-- The previous definitions allow us to calculate excitations into t2g, eg orbitals from 
 +-- either j=1/2 or j=3/2 orbitals. In order to restrict the operator to specific configurations
 +-- we set the restrictions on the transition operator
 +
 +TXASd8x = Clone(TXASx)
 +TXASd9x = Clone(TXASx)
 +TXASd8y = Clone(TXASy)
 +TXASd9y = Clone(TXASy)
 +TXASd8z = Clone(TXASz)
 +TXASd9z = Clone(TXASz)
 +TXASd8x.Restrictions={NF,NB,{"000000 1111111111 0000000000", 9, 9}}
 +TXASd9x.Restrictions={NF,NB,{"000000 1111111111 0000000000",10,10}}
 +TXASd8y.Restrictions={NF,NB,{"000000 1111111111 0000000000", 9, 9}}
 +TXASd9y.Restrictions={NF,NB,{"000000 1111111111 0000000000",10,10}}
 +TXASd8z.Restrictions={NF,NB,{"000000 1111111111 0000000000", 9, 9}}
 +TXASd9z.Restrictions={NF,NB,{"000000 1111111111 0000000000",10,10}}
 +
 +
 +-- We follow the energy definitions as introduced in the group of G.A. Sawatzky (Groningen)
 +-- J. Zaanen, G.A. Sawatzky, and J.W. Allen PRL 55, 418 (1985)
 +-- for parameters of specific materials see
 +-- A.E. Bockquet et al. PRB 55, 1161 (1996)
 +-- After some initial discussion the energies U and Delta refer to the center of a configuration
 +-- The L^10 d^n   configuration has an energy 0
 +-- The L^9  d^n+1 configuration has an energy Delta
 +-- The L^8  d^n+2 configuration has an energy 2*Delta+Udd
 +--
 +-- If we relate this to the onsite energy of the L and d orbitals we find
 +-- 10 eL +  n    ed + n(n-1)     U/2 == 0
 +--  9 eL + (n+1) ed + (n+1)n     U/2 == Delta
 +--  8 eL + (n+2) ed + (n+1)(n+2) U/2 == 2*Delta+U
 +-- 3 equations with 2 unknowns, but with interdependence yield:
 +-- ed = (10*Delta-nd*(19+nd)*U/2)/(10+nd)
 +-- eL = nd*((1+nd)*Udd/2-Delta)/(10+nd)
 +--
 +-- For the final state we/they defined
 +-- The 2p^5 L^10 d^n+1 configuration has an energy 0
 +-- The 2p^5 L^9  d^n+2 configuration has an energy Delta + Udd - Upd
 +-- The 2p^5 L^8  d^n+3 configuration has an energy 2*Delta + 3*Udd - 2*Upd
 +--
 +-- If we relate this to the onsite energy of the p and d orbitals we find
 +-- 6 ep + 10 eL +  n    ed + n(n-1)     Udd/2 + 6 n     Upd == 0
 +-- 6 ep +  9 eL + (n+1) ed + (n+1)n     Udd/2 + 6 (n+1) Upd == Delta
 +-- 6 ep +  8 eL + (n+2) ed + (n+1)(n+2) Udd/2 + 6 (n+2) Upd == 2*Delta+Udd
 +-- 5 ep + 10 eL + (n+1) ed + (n+1)(n)   Udd/2 + 5 (n+1) Upd == 0
 +-- 5 ep +  9 eL + (n+2) ed + (n+2)(n+1) Udd/2 + 5 (n+2) Upd == Delta+Udd-Upd
 +-- 5 ep +  8 eL + (n+3) ed + (n+3)(n+2) Udd/2 + 5 (n+3) Upd == 2*Delta+3*Udd-2*Upd
 +-- 6 equations with 3 unknowns, but with interdependence yield:
 +-- epfinal = (10*Delta + (1+nd)*(nd*Udd/2-(10+nd)*Upd) / (16+nd)
 +-- edfinal = (10*Delta - nd*(31+nd)*Udd/2-90*Upd) / (16+nd)
 +-- eLfinal = ((1+nd)*(nd*Udd/2+6*Upd)-(6+nd)*Delta) / (16+nd)
 +--
 +-- 
 +-- 
 +-- note that ed-ep = Delta - nd * U and not Delta
 +-- note furthermore that ep and ed here are defined for the onsite energy if the system had
 +-- locally nd electrons in the d-shell. In DFT or Hartree Fock the d occupation is in the end not
 +-- nd and thus the onsite energy of the Kohn-Sham orbitals is not equal to ep and ed in model
 +-- calculations.
 +--
 +-- note furthermore that ep and eL actually should be different for most systems. We happily ignore this fact
 +-- 
 +-- We normally take U and Delta as experimentally determined parameters
 +
 +-- number of electrons (formal valence)
 +nd = 8
 +-- parameters from experiment (core level PES)
 +Udd     =  7.3
 +Upd     =  8.5
 +Delta   =  4.7
 +-- parameters obtained from DFT (PRB 85, 165113 (2012))
 +F2dd    = 11.14 
 +F4dd    =  6.87
 +F2pd    =  6.67
 +G1pd    =  4.92
 +G3pd    =  2.80
 +tenDq   =  0.56
 +tenDqL  =  1.44
 +Veg     =  2.06
 +Vt2g    =  1.21
 +zeta_3d =  0.081
 +zeta_2p = 11.51
 +Bz      =  0.000001
 +Hz      =  0.120
 +
 +ed      = (10*Delta-nd*(19+nd)*Udd/2)/(10+nd)
 +eL      = nd*((1+nd)*Udd/2-Delta)/(10+nd)
 +
 +epfinal = (10*Delta + (1+nd)*(nd*Udd/2-(10+nd)*Upd)) / (16+nd)
 +edfinal = (10*Delta - nd*(31+nd)*Udd/2-90*Upd) / (16+nd)
 +eLfinal = ((1+nd)*(nd*Udd/2+6*Upd) - (6+nd)*Delta) / (16+nd)
 +
 +F0dd    = Udd + (F2dd+F4dd) * 2/63
 +F0pd    = Upd + (1/15)*G1pd + (3/70)*G3pd
 +
 +Hamiltonian =  F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d) + Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + ed * OppN_3d + eL * OppN_Ld
 +            
 +XASHamiltonian =  F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d)+ Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + edfinal * OppN_3d + eLfinal * OppN_Ld + epfinal * OppN_2p + zeta_2p * Oppcldots + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3  
 +               
 +-- we now can create the lowest Npsi eigenstates:
 +Npsi=3
 +-- in order to make sure we have a filling of 8 electrons we need to define some restrictions
 +StartRestrictions = {NF, NB, {"000000 1111111111 0000000000",8,8}, {"111111 0000000000 1111111111",16,16}}
 +
 +psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)
 +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
 +
 +-- calculating the spectra is simple and single line once all operators and wave-functions
 +-- are defined.
 +
 +-- we here calculate the spectra for x polarized light and seperate the excitations into the t2g or eg sub-shell.
 +-- We also calculate the spectra separately coming form j=1/2 or j=3/2 orbitals.
 +
 +XASSpectra = CreateSpectra(XASHamiltonian, {TXASx,TXASxeg,TXASxt2g,TXASxj12,TXASxj32,TXASd8x,TXASd9x}, psiList[1], {{"Emin",-10}, {"Emax",20}, {"NE",3000}, {"Gamma",0.1}})
 +
 +XASSpectra.Broaden(0.4, {{-3.7, 0.45}, {-2.2, 0.65}, { 0.0, 0.65}, { 8  , 0.80}, {13.2, 0.80}, {14.0, 1.075}, {16.0, 1.075}})
 +
 +-- in order to plot the spectra we write them to file in ASCII format
 +-- note that the calculated object, i.e. <psi | T^dag 1/(w-H+IG/2) T | psi>
 +-- is complex. Minus the imaginary part is the absorption.
 +XASSpectra.Print({{"file","XASSpec.dat"}})
 +
 +-- in order to plot we use gnuplot. You can use any program, but gnuplot is nice as it can
 +-- be called from within Quanty. We first define the gnuplot script as a string inside Quanty
 +gnuplotInput = [[
 +set autoscale 
 +set xtic auto 
 +set ytic auto 
 +set style line  1 lt 1 lw 1 lc rgb "#000000"
 +set style line  2 lt 1 lw 1 lc rgb "#FF0000"
 +set style line  3 lt 1 lw 1 lc rgb "#008000"
 +
 +set xlabel "E (eV)" font "Times,12"
 +set ylabel "Intensity (arb. units)" font "Times,12"
 +
 +set out "XASpartialSpec.ps"
 +set size 1.0, 1.0
 +set terminal postscript portrait enhanced color  "Times" 12
 +
 +set multiplot layout 5, 1
 +
 +plot "XASSpec.dat"  u 1:(-$3 ) title 'x-total Sz=-1' with lines ls  1,\
 +     "XASSpec.dat"  u 1:(-$5 ) title 'x-eg    Sz=-1' with lines ls  2,\
 +     "XASSpec.dat"  u 1:(-$7 ) title 'x-t2g   Sz=-1' with lines ls  3
 +
 +plot "XASSpec.dat"  u 1:(-$3 ) title 'x total Sz=-1' with lines ls  1,\
 +     "XASSpec.dat"  u 1:(-$9 ) title 'x-j=1/2 Sz=-1' with lines ls  2,\
 +     "XASSpec.dat"  u 1:(-$11) title 'x-j=3/2 Sz=-1' with lines ls  3
 +
 +plot "XASSpec.dat"  u 1:(-$3 ) title 'x total Sz=-1' with lines ls  1,\
 +     "XASSpec.dat"  u 1:(-$13) title 'x-d8    Sz=-1' with lines ls  2,\
 +     "XASSpec.dat"  u 1:(-$15) title 'x-d9    Sz=-1' with lines ls  3
 +
 +]]
 +
 +-- next we save this string to a file
 +file = io.open("XASpartialSpec.gnuplot", "w")
 +file:write(gnuplotInput)
 +file:close()
 +
 +-- and finally call gnuplot to execute the script
 +os.execute("gnuplot XASpartialSpec.gnuplot")
 +-- as I like pdf to view and eps to include in the manuel I transform the format
 +os.execute(" ps2pdf XASpartialSpec.ps; ps2eps XASpartialSpec.ps ;  mv XASpartialSpec.eps temp.eps ; eps2eps temp.eps XASpartialSpec.eps ; rm temp.eps")
 +</code>
 +###
 +
 +The resulting spectra are:
 +| {{ :documentation:tutorials:nio_ligand_field:xaspartialspec.png?nolink |}} |
 +^ $2p$ to $3d$ excitations. Top panel: Split according to excitations into the Ni $3d$ $e_g$ or $t_{2g}$ orbitals. Middle panel: Split according to excitations from the Ni $2p$ $j=1/2$ or $j=3/2$ orbital. Bottom panel: split according to excitations into $d^8$ or $d^9$ configurations. ^
 +
 +###
 +The output of the script is:
 +<file Quanty_Output XAS_partial.out>
 +Start of BlockGroundState. Converge 3 states to an energy with relative variance smaller than  1.490116119384766E-06
 +
 +Start of BlockOperatorPsiSerialRestricted
 +Outer loop   1, Number of Determinants:        45        45 last variance  1.159112471523181E+00
 +Start of BlockOperatorPsiSerialRestricted
 +Start of BlockGroundState. Converge 3 states to an energy with relative variance smaller than  1.490116119384766E-06
 +
 +Start of BlockOperatorPsiSerial
 +Outer loop   1, Number of Determinants:        32       100 last variance  8.359599622001442E+00
 +Start of BlockOperatorPsiSerial
 +Outer loop   2, Number of Determinants:       100       138 last variance  1.315915512292445E+00
 +Start of BlockOperatorPsiSerial
 +  #    <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>
 +  1   -3.503    1.999   12.000   15.095   -0.908   -0.281   -0.305   -1.042   -0.924    2.186    5.990    3.825    6.000    8.175 
 +  2   -3.395    1.999   12.000   15.160   -0.004   -0.002   -0.322   -1.043   -0.925    2.189    5.988    3.823    6.000    8.178 
 +  3   -3.286    1.999   12.000   15.211    0.903    0.278   -0.336   -1.043   -0.925    2.193    5.987    3.820    6.000    8.180 
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Start of LanczosTriDiagonalizeRR
 +Spectra printed to file: XASSpec.dat
 +</file>
 +###
 +
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » NiO ligand field » XAS $L_{2,3}$ partial excitations

Print/export