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.
Navigation
Download area
Workshops
Script versions
External programs
Privacy
Answers
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
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"}}) endcreate_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")