Auger Spectroscopy

asked by Riccardo Piombo (2020/03/19 16:52)

Dear All, Is there any tutorial on how to implement Meitner-Auger process in Quanty? Maybe shown in some past workshop.

Answers

, 2020/03/20 06:45

One can surely do these type of calculations. See Ab initio calculation of the electron capture spectrum of 163Ho: Auger-Meitner decay into continuum states, but there is no comprehensive tutorial available yet.

This will be part of the next workshop in Heidelberg so at least then we will be able to present you with a good tutorial. If you have a specific project in mind and need something earlier feel free to contact us via mail.

Best wishes, Maurits

, 2020/04/13 13:52

By now I am trying to reproduce the d8-DOS of a TM compound in which the TM has an initial electronic structure equal to d10 to see whether my script works correctly. I post here a simplified version of the script in which only a coulomb operator on d shell is taken into account. I'm expecting three states in the singlet spectrum and two in triplet one. But the singlet part gives me five states. Here is the code

NF = 10
NB = 0

nd = 10

IndexDn = {0,2,4,6,8}
IndexUp = {1,3,5,7,9}

A = 7.9   --> +-0.5
B = 0.15
C = 0.58

F0dd = A + 7.0*C/5.0
F2dd = 49.0*B + 7*C
F4dd = 441.0*C/35.0
Udd  = F0dd - (F2dd + F4dd)*2.0/63.0  

OppF0_dd = NewOperator("U", NF, IndexUp, IndexDn, {1,0,0})
OppF0_dd.Name="F^(0) = Coulomb Operator with k=0 in d-shell"
OppF2_dd = NewOperator("U", NF, IndexUp, IndexDn, {0,1,0})
OppF2_dd.Name="F^(2) = Coulomb Operator with k=2 in d-shell"
OppF4_dd = NewOperator("U", NF, IndexUp, IndexDn, {0,0,1})
OppF4_dd.Name="F^(4) = Coulomb Operator with k=4 in d-shell"

print("Coulomb operator for electrons in the d shell")
print(OppF0)
print(OppF2)
print(OppF4)

OppUdd = F0dd*OppF0_dd + F2dd*OppF2_dd + F4dd*OppF4_dd
OppUdd.Name = "dd Coulomb repulsion"

H_r = OppUdd
print("")
print("N-body states energies computation")
StartRestrictions_i = {NF, NB, {"1111111111", nd, nd}}

----------- N-Hamiltonian diagonalization -----------
Npsi_i = 1
psiList_N = Eigensystem(H_r, StartRestrictions_i, Npsi_i)
psiList_N = {psiList_N} 
E_gs_N = psiList_N[1]*(H_r)*psiList_N[1]

----------- change-basis matrix -----------

Orbitals = {"4d"}
YtoK = YtoKMatrix(Orbitals)
KtoY = ConjugateTranspose(YtoK)

----------- (N-2)-Hamiltonian diagonalization -----------

print("")
print("(N-2)-body states energies computation")
StartRestrictions_rem_two_particle = {NF, NB, {"1111111111", nd-2, nd-2}} 

psiList_Nminus2 = Eigensystem(H_r, StartRestrictions_rem_two_particle, 1)
psiList_Nminus2 = {psiList_Nminus2}

E_gs_Nminus2 = psiList_Nminus2[1]*(H_r)*psiList_Nminus2[1]
Ef = E_gs_N - E_gs_Nminus2

----------- Auger TM Spectrum Computation -----------
-- Here you create the operators to put in T_Auger_TM_Singlet and T_Auger_TM_Triplet
dofile("create_Auger_operators.lua") 

E_max = -70
E_min = -150
Gamma =  2.4
Gamma_unb = 0.1
Npoint = 100*(E_max-E_min)/Gamma
Npoint_unb = 100*(E_max-E_min)/Gamma_unb

G_Auger_Singlet = CreateSpectra(H_r, T_Auger_TM_Singlet, psiList_N[1], {{"Emin",E_min}, {"Emax",E_max}, {"NE",Npoint}, {"Gamma",Gamma}} )
G_Auger_Singlet.Shift(Ef)
G_Auger_Singlet = -(1/math.pi)*G_Auger_Singlet

peaks_Auger_Singlet = CreateSpectra(H_r, T_Auger_TM_Singlet, psiList_N[1], {{"Emin", E_min}, {"Emax",E_max}, {"NE",Npoint_unb}, {'Gamma',Gamma_unb}} )
peaks_Auger_Singlet.Shift(Ef)
peaks_Auger_Singlet = -(1/math.pi)*peaks_Auger_Singlet
--------------------------------------------------------------------
G_Auger_Triplet = CreateSpectra(H_r, T_Auger_TM_Triplet, psiList_N[1], {{"Emin",E_min}, {"Emax",E_max}, {"NE",Npoint}, {"Gamma",Gamma}} )
G_Auger_Triplet.Shift(Ef)
G_Auger_Triplet = -(1/math.pi)*G_Auger_Triplet

peaks_Auger_Triplet = CreateSpectra(H_r, T_Auger_TM_Triplet, psiList_N[1], {{"Emin", E_min}, {"Emax",E_max}, {"NE",Npoint_unb}, {'Gamma',Gamma_unb}} )
peaks_Auger_Triplet.Shift(Ef)
peaks_Auger_Triplet = -(1/math.pi)*peaks_Auger_Triplet

if Udd==0 then
  G_Auger_Singlet.Print({{"file", "Auger_Spectrum_TM_Singlet_nonint.dat"}})
  peaks_Auger_Singlet.Print({{"file", "Peaks_Auger_Singlet_nonint.dat"}})
  G_Auger_Triplet.Print({{"file","Auger_Spectrum_TM_Triplet_nonint.dat"}})
  peaks_Auger_Triplet.Print({{"file", "Peaks_Auger_Triplet_nonint.dat"}})
else
  G_Auger_Singlet.Print({{"file", "Auger_Spectrum_TM_Singlet.dat"}})
  peaks_Auger_Singlet.Print({{"file", "Peaks_Auger_Singlet.dat"}})
  G_Auger_Triplet.Print({{"file","Auger_Spectrum_TM_Triplet.dat"}})
  peaks_Auger_Triplet.Print({{"file", "Peaks_Auger_Triplet.dat"}})
end

create_Auger_operators.lua is:

print("")
print("YtoKMatrix without spin")
print("--------------------------------------------")
print(YtoKMatrix({"d"}, {{"addSpin",false}}))

print("")
print("YtoKMatrix conjugate transpose without spin")
print("--------------------------------------------")
print(ConjugateTranspose(YtoKMatrix({"d"}, {{"addSpin",false}})))

print("")
print("YtoKMatrix with spin")
print("--------------------------------------------")
print(YtoK)

print("")
print("YtoKMatrix conjugate transpose with spin")
print("--------------------------------------------")
print(KtoY)

print("")
print("YtoK is the inverse of KtoY and viceversa")
print("--------------------------------------------")
print(Chop(YtoK*KtoY))

--Annihilation Operators with spin up on spherical harmonics basis
op1 = NewOperator("An",NF,1)
op1.Name = "c_-2 up"
op3 = NewOperator("An",NF,3) 
op3.Name = "c_-1 up"
op5 = NewOperator("An",NF,5) 
op5.Name = "c_0 up"
op7 = NewOperator("An",NF,7) 
op7.Name = "c_1 up"
op9 = NewOperator("An",NF,9)  
op9.Name = "c_2 up"

--Annihilation Operators with spin down on spherical harmonics basis
op0 = NewOperator("An",NF,0)
op0.Name = "c_-2 down"
op2 = NewOperator("An",NF,2) 
op2.Name = "c_-1 down"
op4 = NewOperator("An",NF,4) 
op4.Name = "c_0 down"
op6 = NewOperator("An",NF,6) 
op6.Name = "c_+1 down"
op8 = NewOperator("An",NF,8)  
op8.Name = "c_+2 down"

print("")
print("Annihilation operators on a basis |ml> of spherical harmonics")
print("--------------------------------------------")
print("c_-2 up = ") 
print(op1)
print("c_-1 up = ") 
print(op3)
print("c_0 up = ") 
print(op5)
print("c_+1 up = ") 
print(op7)
print("c_+2 up = ") 
print(op9)
print("--------------------------------------------")
print("c_-2 down = ") 
print(op0)
print("c_-1 down = ") 
print(op2)
print("c_0 down = ") 
print(op4)
print("c_+1 down = ") 
print(op6)
print("c_+2 down = ") 
print(op8)


opK1 = Rotate(op1, KtoY) -->  1/sqrt(2)[d_1 + d_9]   =  1/sqrt(2) [ d_{-2 up}   + d_{+2 up} ]  =  d_{x^2-y^2 up}
opK3 = Rotate(op3, KtoY) -->  d_5                    =  d_{0 up}                               =  d_{3z^2-r^2 up}
opK5 = Rotate(op5, KtoY) --> -i/sqrt(2)[ d_3 + d_7 ] = -i/sqrt(2) [ d_{-1 up}   + d_{+1 up} ]  = -d_{yz up}
opK7 = Rotate(op7, KtoY) -->  1/sqrt(2)[ d_3 - d_7 ] =  1/sqrt(2) [ d_{-1 up}   - d_{+1 up} ]  =  d_{xz up}
opK9 = Rotate(op9, KtoY) -->  i/sqrt(2)[ d_9 - d_1 ] =  i/sqrt(2) [ d_{+2 up}   - d_{-2 up} ]  = -d_{xy up}
------------------------------------------------------------------------------------------------------------------------
opK0 = Rotate(op0, KtoY) -->  1/sqrt(2)[d_0 + d_8]   =  1/sqrt(2) [ d_{-2 down} + d_{+2 down}] =  d_{x^2-y^2 down}
opK2 = Rotate(op2, KtoY) -->  d_4                    =  d_{0 down}                             =  d_{3z^2-r^2 down}
opK4 = Rotate(op4, KtoY) --> -i/sqrt(2)[ d_2 + d_6 ] = -i/sqrt(2) [ d_{-1 down} + d_{+1 down}] = -d_{yz down}
opK6 = Rotate(op6, KtoY) -->  1/sqrt(2)[ d_2 - d_6 ] =  1/sqrt(2) [ d_{-1 down} - d_{+1 down}] =  d_{xz down}
opK8 = Rotate(op8, KtoY) -->  i/sqrt(2)[ d_8 - d_0 ] =  i/sqrt(2) [ d_{+2 down} - d_{-2 down}] = -d_{xy down}


print("")
print("Annihilation operator on a basis of cubic harmonics")
print("--------------------------------------------")
print(opK1)
print(opK3)
print(opK5)
print(opK7)
print(opK9)
------------------------------------------
print(opK0)
print(opK2)
print(opK4)
print(opK6)
print(opK8)

print("")
print("two-hole |up,up> triplet states in cubic harmonics basis")
print("--------------------------------------------")

-- up up components
C_yz_upC_xy_up   = (-1)*opK5*(-1)*opK9
C_yz_upC_xy_up.Name = "C_yz_up C_xy_up"
print(C_yz_upC_xy_up)

C_xz_upC_xy_up   = opK7*(-1)*opK9
C_xz_upC_xy_up.Name = "C_xz_up C_xy_up"
print(C_xz_upC_xy_up)

C_x2y2_upC_xy_up = opK1*(-1)*opK9
C_x2y2_upC_xy_up.Name = "C_x2y2_up C_xy_up"
print(C_x2y2_upC_xy_up)

C_3z2_upC_xy_up  = opK3*(-1)*opK9
C_3z2_upC_xy_up.Name = "C_3z2_up C_xy_up"
print(C_3z2_upC_xy_up)
---------------------------------
C_xz_upC_yz_up   = opK7*(-1)*opK5
C_xz_upC_yz_up.Name = "C_xz_up C_yz_up"
print(C_xz_upC_yz_up)

C_x2y2_upC_yz_up = opK1*(-1)*opK5
C_x2y2_upC_yz_up.Name = "C_x2y2_up C_yz_up"
print(C_x2y2_upC_yz_up)

C_3z2_upC_yz_up  = opK3*(-1)*opK5
C_3z2_upC_yz_up.Name = "C_3z2_up C_yz_up"
print(C_3z2_upC_yz_up)
---------------------------------
C_x2y2_upC_xz_up = opK1*opK7
C_x2y2_upC_xz_up.Name = "C_x2y2_up C_xz_up"
print(C_x2y2_upC_xz_up)

C_3z2_upC_xz_up  = opK3*opK7
C_3z2_upC_xz_up.Name = "C_3z2_up C_xz_up"
print(C_3z2_upC_xz_up)
---------------------------------
C_3z2_upC_x2y2_up  = opK3*opK1
C_3z2_upC_x2y2_up.Name = "C_3z2_up C_x2y2_up"
print(C_3z2_upC_x2y2_up)

print("")
print("two-hole |down, down> triplet states in cubic harmonics basis")
print("--------------------------------------------")

-- down down components
C_yz_dnC_xy_dn   = (-1)*opK4*(-1)*opK8
C_yz_dnC_xy_dn.Name = "C_yz_dn C_xy_dn"
print(C_yz_dnC_xy_dn)

C_xz_dnC_xy_dn   = opK6*(-1)*opK8
C_xz_dnC_xy_dn.Name = "C_xz_dn C_xy_dn"
print(C_xz_dnC_xy_dn)

C_x2y2_dnC_xy_dn = opK0*(-1)*opK8
C_x2y2_dnC_xy_dn.Name = "C_x2y2_dn C_xy_dn"
print(C_x2y2_dnC_xy_dn)

C_3z2_dnC_xy_dn  = opK2*(-1)*opK8
C_3z2_dnC_xy_dn.Name = "C_3z2_dn C_xy_dn"
print(C_3z2_dnC_xy_dn)
---------------------------------
C_xz_dnC_yz_dn   = opK6*(-1)*opK4
C_xz_dnC_yz_dn.Name = "C_xz_dn C_yz_dn"
print(C_xz_dnC_yz_dn)

C_x2y2_dnC_yz_dn = opK0*(-1)*opK4
C_x2y2_dnC_yz_dn.Name = "C_x2y2_dn C_yz_dn"
print(C_x2y2_dnC_yz_dn)

C_3z2_dnC_yz_dn  = opK2*(-1)*opK4
C_3z2_dnC_yz_dn.Name = "C_3z2_dn C_yz_dn"
print(C_3z2_dnC_yz_dn)
---------------------------------
C_x2y2_dnC_xz_dn = opK0*opK6
C_x2y2_dnC_xz_dn.Name = "C_x2y2_dn C_xz_dn"
print(C_x2y2_dnC_xz_dn)

C_3z2_dnC_xz_dn  = opK2*opK6
C_3z2_dnC_xz_dn.Name = "C_3z2_dn C_xz_dn"
print(C_3z2_dnC_xz_dn)
---------------------------------
C_3z2_dnC_x2y2_dn  = opK2*opK0
C_3z2_dnC_x2y2_dn.Name = "C_3z2_dn C_x2y2_dn"
print(C_3z2_dnC_x2y2_dn)

print("")
print("two-hole |down, down> triplet states in cubic harmonics basis")
print("--------------------------------------------")

-- |up down> + |down up> components
C_xy_upC_yz_dn__C_xy_dnC_yz_up       = (-1)*opK9*(-1)*opK4 + (-1)*opK8*(-1)*opK5
C_xy_upC_yz_dn__C_xy_dnC_yz_up.Name  = "C_xy_upC_yz_dn + C_xy_dnC_yz_up"
print(C_xy_upC_yz_dn__C_xy_dnC_yz_up)

C_xy_upC_xz_dn__C_xy_dnC_xz_up       = (-1)*opK9*opK6 + (-1)*opK8*opK7
C_xy_upC_xz_dn__C_xy_dnC_xz_up.Name  = "C_xy_upC_xz_dn + C_xy_dnC_xz_up"
print(C_xy_upC_xz_dn__C_xy_dnC_xz_up)

C_xy_upC_x2y2_dn__C_xy_dnC_x2y2_up   = (-1)*opK9*opK0 + (-1)*opK8*opK1
C_xy_upC_x2y2_dn__C_xy_dnC_x2y2_up.Name = "C_xy_upC_x2y2_dn + C_xy_dnC_x2y2_up"
print(C_xy_upC_x2y2_dn__C_xy_dnC_x2y2_up)

C_xy_upC_3z2_dn__C_xy_dnC_3z2_up     = (-1)*opK9*opK2 + (-1)*opK8*opK3
C_xy_upC_3z2_dn__C_xy_dnC_3z2_up.Name = "C_xy_upC_3z2_dn + C_xy_dnC_3z2_up"
print(C_xy_upC_3z2_dn__C_xy_dnC_3z2_up)
------------------------------------------------------
C_yz_upC_xz_dn__C_yz_dnC_xz_up       = (-1)*opK5*opK6 + (-1)*opK4*opK7
C_yz_upC_xz_dn__C_yz_dnC_xz_up.Name  = "C_yz_upC_xz_dn + C_yz_dnC_xz_up" 
print(C_yz_upC_xz_dn__C_yz_dnC_xz_up)

C_yz_upC_x2y2_dn__C_yz_dnC_x2y2_up   = (-1)*opK5*opK0 + (-1)*opK4*opK1
C_yz_upC_x2y2_dn__C_yz_dnC_x2y2_up.Name = "C_yz_upC_x2y2_dn + C_yz_dnC_x2y2_up"
print(C_yz_upC_x2y2_dn__C_yz_dnC_x2y2_up)

C_yz_upC_3z2_dn__C_yz_dnC_3z2_up     = (-1)*opK5*opK2 + (-1)*opK4*opK3
C_yz_upC_3z2_dn__C_yz_dnC_3z2_up.Name = "C_yz_upC_3z2_dn + C_yz_dnC_3z2_up"
print(C_yz_upC_3z2_dn__C_yz_dnC_3z2_up)
------------------------------------------------------
C_xz_upC_x2y2_dn__C_xz_dnC_x2y2_up   = opK7*opK0 + opK6*opK1
C_xz_upC_x2y2_dn__C_xz_dnC_x2y2_up.Name = "C_xz_upC_x2y2_dn + C_xz_dnC_x2y2_up"
print(C_xz_upC_x2y2_dn__C_xz_dnC_x2y2_up)

C_xz_upC_3z2_dn__C_xz_dnC_3z2_up     = opK7*opK2 + opK6*opK3
C_xz_upC_3z2_dn__C_xz_dnC_3z2_up.Name = "C_xz_upC_3z2_dn + C_xz_dnC_3z2_up"
print(C_xz_upC_3z2_dn__C_xz_dnC_3z2_up)
------------------------------------------------------
C_x2y2_upC_3z2_dn__C_x2y2_dnC_3z2_up = opK1*opK2 + opK0*opK3
C_x2y2_upC_3z2_dn__C_x2y2_dnC_3z2_up.Name = "C_x2y2_upC_3z2_dn + C_x2y2_dnC_3z2_up"
print(C_x2y2_upC_3z2_dn__C_x2y2_dnC_3z2_up)

T_Auger_TM_Triplet = {}
T_Auger_TM_Triplet = {C_yz_dnC_xy_dn, C_xz_dnC_xy_dn, C_x2y2_dnC_xy_dn, C_3z2_dnC_xy_dn, C_xz_dnC_yz_dn, C_x2y2_dnC_yz_dn, C_3z2_dnC_yz_dn,
C_x2y2_dnC_xz_dn, C_3z2_dnC_xz_dn, C_3z2_dnC_x2y2_dn, 
C_yz_upC_xy_up, C_xz_upC_xy_up, C_x2y2_upC_xy_up, C_3z2_upC_xy_up, C_xz_upC_yz_up, C_x2y2_upC_yz_up, C_3z2_upC_yz_up, C_x2y2_upC_xz_up,
C_3z2_upC_xz_up, C_3z2_upC_x2y2_up,
C_xy_upC_yz_dn__C_xy_dnC_yz_up, C_xy_upC_xz_dn__C_xy_dnC_xz_up, C_xy_upC_x2y2_dn__C_xy_dnC_x2y2_up, C_xy_upC_3z2_dn__C_xy_dnC_3z2_up, 
C_yz_upC_xz_dn__C_yz_dnC_xz_up, C_yz_upC_x2y2_dn__C_yz_dnC_x2y2_up, C_yz_upC_3z2_dn__C_yz_dnC_3z2_up, C_xz_upC_x2y2_dn__C_xz_dnC_x2y2_up,
C_xz_upC_3z2_dn__C_xz_dnC_3z2_up, C_x2y2_upC_3z2_dn__C_x2y2_dnC_3z2_up}

print("")
print("All two-hole singlet states (|up down> - |down up> ) in cubic harmonics basis")
print("--------------------------------------------")

C_xy_upC_xy_dn   = (-1)*opK9*(-1)*opK8
C_xy_upC_xy_dn.Name = "C_xy_up C_xy_dn"
print(C_xy_upC_xy_dn)

C_xy_upC_yz_dn   = (-1)*opK9*(-1)*opK4
C_xy_upC_yz_dn.Name = "C_xy_up C_yz_dn"
print(C_xy_upC_yz_dn)

C_xy_upC_xz_dn   = (-1)*opK9*opK6
C_xy_upC_xz_dn.Name = "C_xy_up C_xz_dn"
print(C_xy_upC_xz_dn)

C_xy_upC_x2y2_dn = (-1)*opK9*opK0
C_xy_upC_x2y2_dn.Name = "C_xy_up C_x2y2_dn"
print(C_xy_upC_x2y2_dn)

C_xy_upC_3z2_dn  = (-1)*opK9*opK2
C_xy_upC_3z2_dn.Name = "C_xy_up C_3z2_dn"
print(C_xy_upC_3z2_dn)
------------------------------------------------------
C_yz_upC_yz_dn   = (-1)*opK5*(-1)*opK4
C_yz_upC_yz_dn.Name = "C_yz_up C_yz_dn"
print(C_yz_upC_yz_dn)

C_yz_upC_xy_dn   = (-1)*opK5*(-1)*opK8
C_yz_upC_xy_dn.Name = "C_yz_upC_xy_dn"
print(C_yz_upC_xy_dn)

C_yz_upC_xz_dn   = (-1)*opK5*opK6
C_yz_upC_xz_dn.Name = "C_yz_up C_xz_dn"
print(C_yz_upC_xz_dn)

C_yz_upC_x2y2_dn = (-1)*opK5*opK0
C_yz_upC_x2y2_dn.Name = "C_yz_up C_x2y2_dn"
print(C_yz_upC_x2y2_dn)

C_yz_upC_3z2_dn  = (-1)*opK5*opK2
C_yz_upC_3z2_dn.Name = "C_yz_up C_3z2_dn"
print(C_yz_upC_3z2_dn)
------------------------------------------------------
C_xz_upC_xz_dn   = opK7*opK6
C_xz_upC_xz_dn.Name = "C_xz_up C_xz_dn"
print(C_xz_upC_xz_dn)

C_xz_upC_xy_dn   = opK7*(-1)*opK8
C_xz_upC_xy_dn.Name = "C_xz_up C_xy_dn"
print(C_xz_upC_xy_dn)

C_xz_upC_yz_dn   = opK7*(-1)*opK4
C_xz_upC_yz_dn.Name = "C_xz_up C_yz_dn"
print(C_xz_upC_yz_dn)

C_xz_upC_x2y2_dn = opK7*opK0
C_xz_upC_x2y2_dn.Name = "C_xz_up C_x2y2_dn"
print(C_xz_upC_x2y2_dn)

C_xz_upC_3z2_dn  = opK7*opK2
C_xz_upC_3z2_dn.Name = "C_xz_up C_3z2_dn"
print(C_xz_upC_3z2_dn)
------------------------------------------------------
C_x2y2_upC_x2y2_dn = opK1*opK0
C_x2y2_upC_x2y2_dn.Name = "C_x2y2_up C_x2y2_dn"
print(C_x2y2_upC_x2y2_dn)

C_x2y2_upC_xy_dn   = opK1*(-1)*opK8
C_x2ame = "C_x2y2_up C_xy_dn"
print(C_x2y2_upC_xy_dn)

C_x2y2_upC_yz_dn   = opK1*(-1)*opK4
C_x2y2_upC_yz_dn.Name = "C_x2y2_up C_yz_dn"
print(C_x2y2_upC_yz_dn)

C_x2y2_upC_xz_dn   = opK1*opK6
C_x2y2_upC_xz_dn.Name = "C_x2y2_up C_xz_dn"
print(C_x2y2_upC_xz_dn)

C_x2y2_upC_3z2_dn  = opK1*opK2
C_x2y2_upC_3z2_dn.Name = "C_x2y2_up C_3z2_dn"
print(C_x2y2_upC_3z2_dn)
------------------------------------------------------
C_3z2_upC_3z2_dn  = opK3*opK2
C_3z2_upC_3z2_dn.Name = "C_3z2_up C_3z2_dn"
print(C_3z2_upC_3z2_dn)

C_3z2_upC_xy_dn   = opK3*(-1)*opK8
C_3z2_upC_xy_dn.Name = "C_3z2_up C_xy_dn"
print(C_3z2_upC_xy_dn)

C_3z2_upC_yz_dn   = opK3*(-1)*opK4
C_3z2_upC_yz_dn.Name = "C_3z2_up C_yz_dn"
print(C_3z2_upC_yz_dn)

C_3z2_upC_xz_dn   = opK3*opK6
C_3z2_upC_xz_dn.Name = "C_3z2_up C_xz_dn"
print(C_3z2_upC_xz_dn)

C_3z2_upC_x2y2_dn = opK3*opK0
C_3z2_upC_x2y2_dn.Name = "C_3z2_up C_x2y2_dn"
print(C_3z2_upC_x2y2_dn)

T_Auger_TM_Singlet = {}
T_Auger_TM_Singlet = {C_xy_upC_xy_dn, C_xy_upC_yz_dn, C_xy_upC_xz_dn, C_xy_upC_x2y2_dn, C_xy_upC_3z2_dn,
                   C_yz_upC_yz_dn, C_yz_upC_xy_dn, C_yz_upC_xz_dn, C_yz_upC_x2y2_dn, C_yz_upC_3z2_dn,
                   C_xz_upC_xz_dn, C_xz_upC_xy_dn, C_xz_upC_yz_dn, C_xz_upC_x2y2_dn, C_xz_upC_3z2_dn,
                   C_x2y2_upC_x2y2_dn, C_x2y2_upC_xy_dn, C_x2y2_upC_yz_dn, C_x2y2_upC_xz_dn, C_x2y2_upC_3z2_dn,
                   C_3z2_upC_3z2_dn, C_3z2_upC_xy_dn, C_3z2_upC_yz_dn, C_3z2_upC_xz_dn, C_3z2_upC_x2y2_dn}

Then to plot the spectrum use this .py file

from numpy import *
from matplotlib.pyplot import *
from pylab import *
from scipy import integrate
import math, os

# Content of T_Auger_TM_Singlet 
# = {C_xy_upC_xy_dn, C_xy_upC_yz_dn, C_xy_upC_xz_dn, C_xy_upC_x2y2_dn, C_xy_upC_3z2_dn,
#    C_yz_upC_yz_dn, C_yz_upC_xy_dn, C_yz_upC_xz_dn, C_yz_upC_x2y2_dn, C_yz_upC_3z2_dn,
#    C_xz_upC_xz_dn, C_xz_upC_xy_dn, C_xz_upC_yz_dn, C_xz_upC_x2y2_dn, C_xz_upC_3z2_dn,
#    C_x2y2_upC_x2y2_dn, C_x2y2_upC_xy_dn, C_x2y2_upC_yz_dn, C_x2y2_upC_xz_dn, C_x2y2_upC_3z2_dn,
#    C_3z2_upC_3z2_dn, C_3z2_upC_xy_dn, C_3z2_upC_yz_dn, C_3z2_upC_xz_dn, C_3z2_upC_x2y2_dn}


# Content of T_Auger_TM_Triplet
# = {C_yz_dnC_xy_dn, C_xz_dnC_xy_dn, C_x2y2_dnC_xy_dn, C_3z2_dnC_xy_dn, C_xz_dnC_yz_dn, C_x2y2_dnC_yz_dn, C_3z2_dnC_yz_dn,
#    C_x2y2_dnC_xz_dn, C_3z2_dnC_xz_dn, C_3z2_dnC_x2y2_dn, 
#    C_yz_upC_xy_up, C_xz_upC_xy_up, C_x2y2_upC_xy_up, C_3z2_upC_xy_up, C_xz_upC_yz_up, C_x2y2_upC_yz_up, C_3z2_upC_yz_up, C_x2y2_upC_xz_up,
#    C_3z2_upC_xz_up, C_3z2_upC_x2y2_up,
#    C_xy_upC_yz_dn__C_xy_dnC_yz_up, C_xy_upC_xz_dn__C_xy_dnC_xz_up, C_xy_upC_x2y2_dn__C_xy_dnC_x2y2_up, C_xy_upC_3z2_dn__C_xy_dnC_3z2_up, 
#    C_yz_upC_xz_dn__C_yz_dnC_xz_up, C_yz_upC_x2y2_dn__C_yz_dnC_x2y2_up, C_yz_upC_3z2_dn__C_yz_dnC_3z2_up, C_xz_upC_x2y2_dn__C_xz_dnC_x2y2_up,
#    C_xz_upC_3z2_dn__C_xz_dnC_3z2_up, C_x2y2_upC_3z2_dn__C_x2y2_dnC_3z2_up}

################## d sector ##################
#----------------------------------------------------------------
# what is yet to be done is to sum up all the columns

# The program works as follow:
# reduce an RxC matrix to an Rx2 matrix where the second column is  
# the sum of the (C-1) elements in columns 2, .... C of the first matrix

print("Loading Auger Singlet")
Auger_Singlet = loadtxt("Auger_Spectrum_TM_Singlet.dat", skiprows=1)
sumbyrow1 = column_stack((Auger_Singlet[:,0],sum(Auger_Singlet[:,1:],axis=1)))
print("Done!")

print("Loading Auger Triplet")
Auger_Triplet = loadtxt("Auger_Spectrum_TM_Triplet.dat", skiprows=1)
sumbyrow2 = column_stack((Auger_Triplet[:,0],sum(Auger_Triplet[:,1:],axis=1)))
print("Done!")

w = Auger_Singlet[:,0]

print("Loading Auger Singlet Peaks")
Auger_Singlet_Peaks = loadtxt("Peaks_Auger_Singlet.dat", skiprows=1)
sumbyrow3 = column_stack((Auger_Singlet_Peaks[:,0],sum(Auger_Singlet_Peaks[:,1:],axis=1)))
print("Done!")

print("Loading Auger Triplet Peaks")
Auger_Triplet_Peaks = loadtxt("Peaks_Auger_Triplet.dat", skiprows=1)
sumbyrow4 = column_stack((Auger_Triplet_Peaks[:,0],sum(Auger_Triplet_Peaks[:,1:],axis=1)))
print("Done!")

wp = Auger_Singlet_Peaks[:,0]

Total_Auger = sumbyrow1[:,1] + sumbyrow2[:,1]
Total_Auger_Peaks = sumbyrow3[:,1] + sumbyrow4[:,1]

max1 = amax(Total_Auger)
max2 = amax(Total_Auger_Peaks)

#print("Loading Auger experimental data")
#Auger_exp = loadtxt("AgF_SHIFT_Auger.dat")
#w_exp = Auger_exp[:,0]
#experiment = Auger_exp[:,1]
#max_exp = amax(experiment)
#print("Done!")

# each submyrow contains 2 columns: in the first column there are the frequencies of the 
# two-body GF while in the second one there is the total GF

Gamma = 0.1

figure()
xlim(20,-20)
plot(sumbyrow1[:,0], sumbyrow1[:,1], label = 'Auger Singlet', color = 'b')
plot(sumbyrow3[:,0], sumbyrow3[:,1]*Gamma/10, label = 'Auger Singlet Peaks', color = 'b')
plot(sumbyrow2[:,0], sumbyrow2[:,1]/3., label = 'Auger Triplet', color = 'g', linestyle = 'dashed')
plot(sumbyrow4[:,0], sumbyrow4[:,1]*Gamma/10, label = 'Auger Triplet Peaks', color = 'g', linestyle = 'dashed')
legend()

figure()
xlim(20,-20)
plot(w,Total_Auger/max1, label = 'Total Auger')
plot(wp,Total_Auger_Peaks*Gamma/max2, label = 'Total Auger Peaks')
#plot(w_exp, experiment/max_exp+0.5, label="experimental MNN Auger")
You could leave a comment if you were logged in.
Print/export