In order to shorten the examples, make them better readable and focus more on the new part introduced in the calculation we split the input file into several parts which are saved as separate files. The input file attached here defines all operators and calculates the ground-state. The tutorials use these to calculate spectra.
-- The example on Fluorescence yield, RIXS and non-resonant IXS all need a very similar
-- starting point. We need to define the Ground-state and some operators on this state.
-- We need to define Hamiltonians with either a 2p or 3s core hole and interactions between
-- these.
-- Here we create a file that is included in tutorials later
Verbosity(0)
NF=28
NB=0
IndexDn_2p={ 0, 2, 4}
IndexUp_2p={ 1, 3, 5}
IndexDn_3s={ 6}
IndexUp_3s={ 7}
IndexDn_3d={ 8,10,12,14,16}
IndexUp_3d={ 9,11,13,15,17}
IndexDn_Ld={18,20,22,24,26}
IndexUp_Ld={19,21,23,25,27}
-- 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_3s = NewOperator("Number", NF, IndexUp_3s, IndexUp_3s, {1})
OppNDn_3s = NewOperator("Number", NF, IndexDn_3s, IndexDn_3s, {1})
OppN_3s = OppNUp_3s + OppNDn_3s
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 2p to 3d
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})
-- core valence interaction 3s to 3d
OppUsdF0 = NewOperator("U", NF, IndexUp_3s, IndexDn_3s, IndexUp_3d, IndexDn_3d, {1}, {0})
OppUsdG2 = NewOperator("U", NF, IndexUp_3s, IndexDn_3s, IndexUp_3d, IndexDn_3d, {0}, {1})
-- next we define the dipole operator. The dipole operator is given as epsilon.r
-- with epsilon the polarization vector of the light and r the unit position vector
-- We can expand the position vector on (renormalized) spherical harmonics and use
-- the crystal-field operator to create the dipole operator.
-- we both define the dipole operator between the 2p and 3d shell as well as the dipole
-- operator between the 3s and 2p shell
-- x polarized light is defined as x = Cos[phi]Sin[theta] = sqrt(1/2) ( C_1^{(-1)} - C_1^{(1)})
Akm = {{1,-1,sqrt(1/2)},{1, 1,-sqrt(1/2)}}
T2p3dx = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
T3s2px = NewOperator("CF", NF, IndexUp_2p, IndexDn_2p, IndexUp_3s, IndexDn_3s, Akm)
-- y polarized light is defined as y = Sin[phi]Sin[theta] = sqrt(1/2) I ( C_1^{(-1)} + C_1^{(1)})
Akm = {{1,-1,sqrt(1/2)*I},{1, 1,sqrt(1/2)*I}}
T2p3dy = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
T3s2py = NewOperator("CF", NF, IndexUp_2p, IndexDn_2p, IndexUp_3s, IndexDn_3s, Akm)
-- z polarized light is defined as z = Cos[theta] = C_1^{(0)}
Akm = {{1,0,1}}
T2p3dz = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
T3s2pz = NewOperator("CF", NF, IndexUp_2p, IndexDn_2p, IndexUp_3s, IndexDn_3s, Akm)
-- besides linear polarized light one can define circular polarized light as the sum of
-- x and y polarizations with complex prefactors
T2p3dr = sqrt(1/2)*(T2p3dx - I * T2p3dy)
T2p3dl =-sqrt(1/2)*(T2p3dx + I * T2p3dy)
-- we can remove zero's from the dipole operator by chopping it.
T2p3dr.Chop()
T2p3dl.Chop()
-- the 3d to 2p dipole transition is the conjugate transpose of the 2p to 3d dipole transition
T3d2px = ConjugateTranspose(T2p3dx)
T3d2py = ConjugateTranspose(T2p3dy)
T3d2pz = ConjugateTranspose(T2p3dz)
T3d2pl = ConjugateTranspose(T2p3dl)
T3d2pr = ConjugateTranspose(T2p3dr)