LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-06 22:24:03 +01:00
Code Releases Activity

T-matrix self-energy

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-10-06 22:35:36 +02:00
parent 6ca04d4b44
commit a5047e9b1a
11 changed files with 429 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF MOM
T F F
# MP2 MP3 MP2-F12
F F F
T F F
# CCD CCSD CCSD(T)
F F F
# CIS TDHF ppRPA ADC

2
input/options
View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.00001 F 1
# CIS/TDHF/BSE: singlet triplet
T T
T F
# GF: maxSCF thresh DIIS n_diis renormalization
64 0.00001 T 10 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta

2
src/QuAcK/QuAcK.f90
View File

@ -408,7 +408,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA

77
src/QuAcK/excitation_density_Tmatrix.f90 Normal file
View File

@ -0,0 +1,77 @@
subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,X2,Y2,rho1,rho2)
! Compute excitation densities for T-matrix self-energy
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nOO,nVV
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(out) :: X1(nVV,nVV)
double precision,intent(out) :: Y1(nOO,nVV)
double precision,intent(out) :: X2(nVV,nOO)
double precision,intent(out) :: Y2(nOO,nOO)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
integer :: p
integer :: ab,cd,ij,kl
! Output variables
double precision,intent(out) :: rho1(nBas,nBas,nVV)
double precision,intent(out) :: rho2(nBas,nBas,nOO)
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVV
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
enddo
enddo
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
enddo
enddo
enddo
enddo
do a=nO+1,nBas-nR
do ij=1,nOO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
enddo
enddo
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,a,ij) = rho1(p,a,ij) + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
enddo
enddo
enddo
enddo
enddo
end subroutine excitation_density_Tmatrix

88
src/QuAcK/linear_response_ph.f90 Normal file
View File

@ -0,0 +1,88 @@
subroutine linear_response_ph(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_phRPA,Omega,XpY)
! Compute the p-h channel of the linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS)
! Local variables
double precision,external :: trace_matrix
double precision,allocatable :: A(:,:),B(:,:),M(:,:),w(:)
! Output variables
double precision,intent(out) :: Ec_phRPA
double precision,intent(out) :: Omega(nS),XpY(nS,nS)
! Memory allocation
allocate(A(nS,nS),B(nS,nS),M(2*nS,2*nS),w(2*nS))
! Build A and B matrices
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A)
! Tamm-Dancoff approximation
B(:,:) = 0d0
if(.not. TDA) then
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B)
endif
!------------------------------------------------------------------------
! Solve the p-h eigenproblem
!------------------------------------------------------------------------
!
! | +A +B | | X Y | | w 0 | | X Y |
! | | | | = | | | |
! | -B -A | | Y X | | 0 -w | | Y X |
!
! Diagonal blocks
M(1:nS,1:nS) = +A(1:nS,1:nS)
M(nS+1:2*nS,nS+1:2*nS) = -A(1:nS,1:nS)
! Off-diagonal blocks
M(1:nS,nS+1:2*nS) = -B(1:nS,1:nS)
M(nS+1:2*nS,1:nS) = +B(1:nS,1:nS)
! Diagonalize the p-h matrix
call diagonalize_matrix(2*nS,M(:,:),w(:))
Omega(1:nS) = w(nS+1:2*nS)
! Build X+Y
XpY(1:nS,1:nS) = M(nS+1:2*nS,1:nS) + M(nS+1:2*nS,nS+1:2*nS)
call DA(nS,1d0/sqrt(Omega),XpY)
! print*,'X+Y'
! call matout(nS,nS,XpY)
! print*,'RPA excitations'
call matout(2*nS,1,w(:))
! Compute the RPA correlation energy
Ec_phRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
print*,'Ec_phRPA = ',Ec_phRPA
end subroutine linear_response_ph

63
src/QuAcK/linear_response_pp.f90
View File

@ -1,4 +1,5 @@
subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_ppRPA)
subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA)
! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013)
@ -10,42 +11,35 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E
logical,intent(in) :: dRPA
logical,intent(in) :: TDA
logical,intent(in) :: BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: nOO
integer :: nVV
double precision :: trace_matrix
double precision,allocatable :: B(:,:)
double precision,allocatable :: C(:,:)
double precision,allocatable :: D(:,:)
double precision,allocatable :: M(:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: w(:)
double precision,allocatable :: w1(:)
double precision,allocatable :: w2(:)
double precision,allocatable :: X1(:,:)
double precision,allocatable :: Y1(:,:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: Y2(:,:)
double precision,allocatable :: Omega(:)
! Output variables
double precision,intent(out) :: Omega1(nVV)
double precision,intent(out) :: X1(nVV,nVV)
double precision,intent(out) :: Y1(nOO,nVV)
double precision,intent(out) :: Omega2(nOO)
double precision,intent(out) :: X2(nVV,nOO)
double precision,intent(out) :: Y2(nOO,nOO)
double precision,intent(out) :: Ec_ppRPA
! Useful quantities
nOO = nO*(nO-1)/2
nVV = nV*(nV-1)/2
! Memory allocation
allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),w(nOO+nVV), &
w1(nVV),w2(nOO),X1(nVV,nVV),Y1(nVV,nOO),X2(nOO,nVV),Y2(nOO,nOO))
allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV))
! Build B, C and D matrices for the pp channel
@ -78,32 +72,31 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E
! Diagonalize the p-h matrix
Z(:,:) = M(:,:)
call diagonalize_matrix(nOO+nVV,Z(:,:),w(:))
call diagonalize_matrix(nOO+nVV,Z(:,:),Omega(:))
write(*,*) 'pp-RPA excitation energies'
call matout(nOO+nVV,1,w(:))
write(*,*)
! write(*,*) 'pp-RPA excitation energies'
! call matout(nOO+nVV,1,Omega(:))
! write(*,*)
! Split the various quantities in p-p and h-h parts
w1(:) = w(nOO+1:nOO+nVV)
w2(:) = w(1:nOO)
Omega1(:) = Omega(nOO+1:nOO+nVV)
Omega2(:) = Omega(1:nOO)
X1(:,:) = Z(nOO+1:nOO+nVV, 1:nVV )
Y1(:,:) = Z(nOO+1:nOO+nVV,nVV+1:nOO+nVV)
X2(:,:) = Z( 1:nOO , 1:nVV )
Y2(:,:) = Z( 1:nOO ,nVV+1:nOO+nVV)
X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV)
Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO )
Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
if(minval(w1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(w2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!')
if(minval(Omega1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(Omega2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!')
! Compute the RPA correlation energy
Ec_ppRPA = 0.5d0*( sum(w1(:)) - sum(w2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
Ec_ppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
print*,'Ec(pp-RPA) = ',Ec_ppRPA
print*,'Ec(pp-RPA) = ',0.5d0*( sum(abs(w(:))) - trace_matrix(nVV*nOO,M(:,:)))
print*,'Ec(pp-RPA) = ',+sum(w1(:)) - trace_matrix(nVV,C(:,:))
print*,'Ec(pp-RPA) = ',-sum(w2(:)) - trace_matrix(nOO,D(:,:))
print*,'Ec(pp-RPA) = ',+sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
print*,'Ec(pp-RPA) = ',-sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
end subroutine linear_response_pp

40
src/QuAcK/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -14,7 +14,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: e(nBas)
@ -26,10 +25,15 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
logical :: TDA
logical :: BSE
integer :: ispin
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
integer :: nOO
integer :: nVV
double precision,allocatable :: Omega1(:,:)
double precision,allocatable :: X1(:,:,:)
double precision,allocatable :: Y1(:,:,:)
double precision,allocatable :: Omega2(:,:)
double precision,allocatable :: X2(:,:,:)
double precision,allocatable :: Y2(:,:,:)
double precision :: rho
double precision :: Ec_ppRPA(nspin)
! Hello world
@ -40,6 +44,11 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
write(*,*)'****************************************'
write(*,*)
! Useful quantities
nOO = nO*(nO-1)/2
nVV = nV*(nV-1)/2
! Initialization
Ec_ppRPA(:) = 0d0
@ -58,7 +67,8 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
! Singlet manifold
@ -66,9 +76,12 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
ispin = 1
call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, &
Ec_ppRPA)
! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin))
call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
endif
@ -78,9 +91,12 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF
ispin = 2
call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, &
Ec_ppRPA)
! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin))
call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
endif

4
src/QuAcK/print_excitation.f90
View File

@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
! Input variables
character*5,intent(in) :: method
character*12,intent(in) :: method
integer,intent(in) :: ispin,nS
double precision,intent(in) :: Omega(nS)
@ -22,7 +22,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
write(*,*)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A4,A14,A7,A9,A25)')'|',method,' calculation: ',spin_manifold,' manifold',' |'
write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold',' |'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'

65
src/QuAcK/renormalization_factor_Tmatrix.f90 Normal file
View File

@ -0,0 +1,65 @@
subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,Z)
! Compute renormalization factor of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) - 2d0*rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
eps = e(p) + e(a) - Omega2(kl)
Z(p) = Z(p) - 2d0*rho2(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
enddo
! Compute renormalization factor from derivative of SigT
Z(:) = 1d0/(1d0 - Z(:))
end subroutine renormalization_factor_Tmatrix

71
src/QuAcK/self_energy_Tmatrix.f90 Normal file
View File

@ -0,0 +1,71 @@
subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,SigT)
! Compute the correlation part of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,j,k,l,a,b,c,d,p,q,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: SigT(nBas,nBas)
! Initialize
SigT = 0d0
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
SigT(p,q) = SigT(p,q) + 2d0*rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
eps = e(p) + e(a) - Omega2(kl)
SigT(p,q) = SigT(p,q) + 2d0*rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
enddo
end subroutine self_energy_Tmatrix

67
src/QuAcK/self_energy_Tmatrix_diag.f90 Normal file
View File

@ -0,0 +1,67 @@
subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,SigT)
! Compute diagonal of the correlation part of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: SigT(nBas)
! Initialize
SigT = 0d0
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + 2d0*rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
eps = e(p) + e(a) - Omega2(kl)
SigT(p) = SigT(p) + 2d0*rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
end subroutine self_energy_Tmatrix_diag