PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 15:14:55 +02:00
Code Releases Activity

trying a sanity check

Browse Source
This commit is contained in:
Antoine Marie 2025-04-21 14:56:27 +02:00
parent c3d222ee26
commit 921902962b
4 changed files with 434 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/GT/GG0T0pp.f90
View File

@ -128,7 +128,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
if(regularize) call GTpp_regularization(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z,ERI)
!----------------------------------------------
! Solve the quasi-particle equation

101
src/GT/GGTpp_self_energy_diag.f90
View File

@ -1,4 +1,4 @@
subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z,ERI)
! Compute diagonal of the correlation part of the T-matrix self-energy
@ -20,11 +20,12 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b,p,cd,kl
double precision :: num,eps
integer :: i,j,k,a,b,c,p,m,cd,kl
double precision :: num,eps,dem1,dem2
! Output variables
@ -72,6 +73,100 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh
end do
end do
!-----------------------------------------------!
! Testing another way to compute GT self-energy !
!-----------------------------------------------!
! do p=nC+1,nBas-nR
! do i=nC+1,nO
! do j=nC+1,nO
! do a=nO+1,nBas-nR
! eps = e(p) + e(a) - e(i) - e(j)
! num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2
! Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
! Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
! end do
! do a=nO+1,nBas-nR
! do m=1,nVV
! num = - ERI(p,a,i,j) * rho1(p,a,m) * rho1(i,j,m)
! dem1 = e(p) + e(a) - e(i) - e(j)
! dem2 = Om1(m) - e(i) - e(j)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2
! end do
! do m=1,nOO
! num = - ERI(p,a,i,j) * rho2(p,a,m) * rho2(i,j,m)
! dem1 = e(p) + e(a) - e(i) - e(j)
! dem2 = e(p) + e(a) - Om2(m)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2
! end do
! end do
! do k=nC+1,nO
! do m=1,nVV
! num = - ERI(p,i,j,k) * rho1(p,i,m) * rho1(j,k,m)
! dem1 = e(p) + e(i) - Om1(m)
! dem2 = Om1(m) - e(j) - e(k)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2
! end do
! end do
! end do
! end do
! end do
! do p=nC+1,nBas-nR
! do a=nO+1,nBas-nR
! do b=nO+1,nBas-nR
! do i=nC+1,nO
! eps = e(p) + e(i) - e(a) - e(b)
! num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2
! Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
! Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
! end do
! do i=nC+1,nO
! do m=1,nVV
! num = ERI(p,i,a,b) * rho1(p,i,m) * rho1(a,b,m)
! dem1 = e(p) + e(i) - e(a) - e(b)
! dem2 = e(p) + e(i) - Om1(m)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2
! end do
! do m=1,nOO
! num = ERI(p,i,a,b) * rho2(p,i,m) * rho2(a,b,m)
! dem1 = e(p) + e(i) - e(a) - e(b)
! dem2 = Om2(m) - e(a) - e(b)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2
! end do
! end do
! do c=nO+1,nBas-nR
! do m=1,nOO
! num = ERI(p,a,b,c) * rho2(p,a,m) * rho2(b,c,m)
! dem1 = e(p) + e(a) - Om2(m)
! dem2 = Om2(m) - e(b) - e(c)
! Sig(p) = Sig(p) + num/dem1/dem2
! Z(p) = Z(p) - num/dem1/dem1/dem2
! end do
! end do
! end do
! end do
! end do
!-------------------------------------!
! Galitskii-Migdal correlation energy !
!-------------------------------------!

264
src/Parquet/G_Parquet_self_energy.f90
View File

@ -89,127 +89,151 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
! eh part of the self-energy !
!-----------------------------!
call wall_time(start_t)
! call wall_time(start_t)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
!$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_rho,eh_Om,SigC,Z)
!$OMP DO COLLAPSE(2)
do p=nC+1,nOrb-nR
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
! !$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_rho,eh_Om,SigC,Z)
! !$OMP DO COLLAPSE(2)
! do p=nC+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nOrb-nR
! do i=nC+1,nO
! do a=nO+1,nOrb-nR
do n=1,nS
!3h2p
do j=nC+1,nO
num = - ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) &
+ ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n)
dem1 = eQP(a) - eQP(i) - eh_Om(n)
dem2 = eQP(p) - eQP(j) + eh_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! do n=1,nS
! !3h2p
! do j=nC+1,nO
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = - ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) &
! + ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n)
! dem1 = eQP(p) - eQP(j) + eh_Om(n)
! dem2 = eQP(p) - eQP(j) + eQP(a) - eQP(i)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) &
! - num * (reg1/dem1/dem1) * (reg2/dem2)
num = + ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) &
- ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n)
! ! num = - 0d0*ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) &
! ! + ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n)
! ! dem1 = eQP(a) - eQP(i) - eh_Om(n)
! ! dem2 = eQP(p) - eQP(j) + eh_Om(n)
! ! reg1 = (1d0 - 0d0*exp(- 2d0 * eta * dem1 * dem1))
! ! reg2 = (1d0 - 0d0*exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = - ERI(p,i,j,a) * eh_rho(p,j,nS+n) * eh_rho(a,i,nS+n) &
+ ERI(p,i,a,j) * eh_rho(i,a,n) * eh_rho(j,p,n)
! ! num = + ERI(p,a,j,i) * eh_rho(p,j,nS+n) * eh_rho(i,a,nS+n) &
! ! - ERI(p,a,i,j) * eh_rho(a,i,n) * eh_rho(j,p,n)
dem1 = eQP(a) - eQP(i) + eh_Om(n)
dem2 = eQP(p) - eQP(j) + eh_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! ! dem1 = eQP(a) - eQP(i) - eh_Om(n)
! ! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! ! reg1 = (1d0 - 0d0*exp(- 2d0 * eta * dem1 * dem1))
! ! reg2 = (1d0 - 0d0*exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = - ERI(p,i,j,a) * eh_rho(p,j,nS+n) * eh_rho(a,i,nS+n) &
! + ERI(p,i,a,j) * eh_rho(i,a,n) * eh_rho(j,p,n)
! dem1 = eQP(a) - eQP(i) + eh_Om(n)
! dem2 = eQP(p) - eQP(j) + eh_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = - ERI(p,a,j,i) * eh_rho(p,j,n) * eh_rho(i,a,n) &
+ ERI(p,a,i,j) * eh_rho(a,i,nS+n) * eh_rho(j,p,nS+n)
! num = - ERI(p,a,j,i) * eh_rho(p,j,n) * eh_rho(i,a,n) &
! + ERI(p,a,i,j) * eh_rho(a,i,nS+n) * eh_rho(j,p,nS+n)
dem1 = eQP(a) - eQP(i) + eh_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! dem1 = eQP(a) - eQP(i) + eh_Om(n)
! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = - ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) &
+ ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n)
! end do ! j
! !3p2h
! do b=nO+1,nOrb-nR
! num = - ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) &
! + ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n)
dem1 = eQP(a) - eQP(i) - eh_Om(n)
dem2 = eQP(p) - eQP(b) - eh_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! dem1 = eQP(p) - eQP(b) - eh_Om(n)
! dem2 = eQP(p) - eQP(b) - eQP(a) + eQP(i)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) &
! - num * (reg1/dem1/dem1) * (reg2/dem2)
num = + ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) &
- ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n)
! ! num = - ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) &
! ! + ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n)
dem1 = eQP(a) - eQP(i) - eh_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! ! dem1 = eQP(a) - eQP(i) - eh_Om(n)
! ! dem2 = eQP(p) - eQP(b) - eh_Om(n)
! ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = - ERI(p,a,b,i) * eh_rho(p,b,n) * eh_rho(i,a,n) &
+ ERI(p,a,i,b) * eh_rho(a,i,nS+n) * eh_rho(b,p,nS+n)
! ! num = + ERI(p,i,b,a) * eh_rho(p,b,n) * eh_rho(a,i,n) &
! ! - ERI(p,i,a,b) * eh_rho(i,a,nS+n) * eh_rho(b,p,nS+n)
dem1 = eQP(a) - eQP(i) + eh_Om(n)
dem2 = eQP(p) - eQP(b) - eh_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! ! dem1 = eQP(a) - eQP(i) - eh_Om(n)
! ! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! ! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! ! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = - ERI(p,i,b,a) * eh_rho(p,b,nS+n) * eh_rho(a,i,nS+n) &
+ ERI(p,i,a,b) * eh_rho(i,a,n) * eh_rho(b,p,n)
! num = - ERI(p,a,b,i) * eh_rho(p,b,n) * eh_rho(i,a,n) &
! + ERI(p,a,i,b) * eh_rho(a,i,nS+n) * eh_rho(b,p,nS+n)
dem1 = eQP(a) - eQP(i) + eh_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! dem1 = eQP(a) - eQP(i) + eh_Om(n)
! dem2 = eQP(p) - eQP(b) - eh_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! b
! num = - ERI(p,i,b,a) * eh_rho(p,b,nS+n) * eh_rho(a,i,nS+n) &
! + ERI(p,i,a,b) * eh_rho(i,a,n) * eh_rho(b,p,n)
! dem1 = eQP(a) - eQP(i) + eh_Om(n)
! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! end do ! b
end do ! n
! end do ! n
end do ! a
end do ! i
! end do ! a
! end do ! i
end do ! p
!$OMP END DO
!$OMP END PARALLEL
! end do ! p
! !$OMP END DO
! !$OMP END PARALLEL
call wall_time(end_t)
t = end_t - start_t
! call wall_time(end_t)
! t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds'
write(*,*)
! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds'
! write(*,*)
!-----------------------------!
! pp part of the self-energy !
@ -257,22 +281,32 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
do c=nO+1,nOrb-nR
num = - ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
dem1 = hh_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - hh_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
dem1 = hh_Om(n) - eQP(i) - eQP(j)
dem1 = eQP(p) + eQP(c) - hh_Om(n)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) &
- num * (reg1/dem1/dem1) * (reg2/dem2)
! num = - ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
! dem1 = hh_Om(n) - eQP(i) - eQP(j)
! dem2 = eQP(p) + eQP(c) - hh_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
! dem1 = hh_Om(n) - eQP(i) - eQP(j)
! dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! c
end do ! n
@ -322,22 +356,32 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
do k=nC+1,nO
num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
dem1 = ee_Om(n) - eQP(a) - eQP(b)
dem1 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - ee_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2) &
- num * (reg1/dem1/dem1) * (reg2/dem2)
! num = ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
! dem1 = ee_Om(n) - eQP(a) - eQP(b)
! dem2 = eQP(p) + eQP(k) - ee_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
num = - ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
dem1 = ee_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = - ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
! dem1 = ee_Om(n) - eQP(a) - eQP(b)
! dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! c
end do ! n

362
src/Parquet/R_Parquet_self_energy.f90
View File

@ -89,249 +89,249 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
!-------------------------------------!
! singlet eh part of the self-energy !
!-------------------------------------!
! call wall_time(start_t)
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
! !$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_sing_rho,eh_sing_Om,SigC,Z)
! !$OMP DO COLLAPSE(2)
! do p=nC+1,nOrb-nR
call wall_time(start_t)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
!$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_sing_rho,eh_sing_Om,SigC,Z)
!$OMP DO COLLAPSE(2)
do p=nC+1,nOrb-nR
! do i=nC+1,nO
! do a=nO+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nOrb-nR
! do n=1,nS
! !3h2p
! do j=nC+1,nO
! num = ( - ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* &
! eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n)
do n=1,nS
!3h2p
do j=nC+1,nO
num = ( - ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* &
eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n)
! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
! dem2 = eQP(p) - eQP(j) + eh_sing_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
dem2 = eQP(p) - eQP(j) + eh_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (ERI(p,a,j,i) - 0.5d0*ERI(p,a,i,j))* &
! eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n)
num = (ERI(p,a,j,i) - 0.5d0*ERI(p,a,i,j))* &
eh_sing_rho(a,i,n) * eh_sing_rho(j,p,n)
! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (- ERI(p,i,j,a) + 0.5d0*ERI(p,i,a,j)) * &
! eh_sing_rho(j,p,n) * eh_sing_rho(i,a,n)
num = (- ERI(p,i,j,a) + 0.5d0*ERI(p,i,a,j)) * &
eh_sing_rho(j,p,n) * eh_sing_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
! dem2 = eQP(p) - eQP(j) + eh_sing_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
dem2 = eQP(p) - eQP(j) + eh_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (- ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* &
! eh_sing_rho(p,j,n) * eh_sing_rho(i,a,n)
num = (- ERI(p,a,j,i) + 0.5d0*ERI(p,a,i,j))* &
eh_sing_rho(p,j,n) * eh_sing_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! end do ! j
! !3p2h
! do b=nO+1,nOrb-nR
! num = (- ERI(p,a,b,i) + 0.5d0*ERI(p,a,i,b)) * &
! eh_sing_rho(p,b,n) * eh_sing_rho(i,a,n)
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = (- ERI(p,a,b,i) + 0.5d0*ERI(p,a,i,b)) * &
eh_sing_rho(p,b,n) * eh_sing_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
! dem2 = eQP(p) - eQP(b) - eh_sing_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
dem2 = eQP(p) - eQP(b) - eh_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * &
! eh_sing_rho(b,p,n) * eh_sing_rho(i,a,n)
num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * &
eh_sing_rho(b,p,n) * eh_sing_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_sing_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * &
! eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n)
num = (- ERI(p,i,b,a) + 0.5d0*ERI(p,i,a,b)) * &
eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n)
! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
! dem2 = eQP(p) - eQP(b) - eh_sing_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
dem2 = eQP(p) - eQP(b) - eh_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = (ERI(p,i,b,a) - 0.5d0*ERI(p,i,a,b)) * &
! eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n)
num = (ERI(p,i,b,a) - 0.5d0*ERI(p,i,a,b)) * &
eh_sing_rho(p,b,n) * eh_sing_rho(a,i,n)
! dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_sing_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! end do ! b
end do ! b
! end do ! n
end do ! n
! end do ! a
! end do ! i
end do ! a
end do ! i
! end do ! p
! !$OMP END DO
! !$OMP END PARALLEL
! call wall_time(end_t)
! t = end_t - start_t
end do ! p
!$OMP END DO
!$OMP END PARALLEL
call wall_time(end_t)
t = end_t - start_t
! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh self-energy =',t,' seconds'
! write(*,*)
! !-------------------------------------!
! ! triplet eh part of the self-energy !
! !-------------------------------------!
! call wall_time(start_t)
! !$OMP PARALLEL DEFAULT(NONE) &
! !$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
! !$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_trip_rho,eh_trip_Om,SigC,Z)
! !$OMP DO COLLAPSE(2)
! do p=nC+1,nOrb-nR
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for singlet eh self-energy =',t,' seconds'
write(*,*)
!-------------------------------------!
! triplet eh part of the self-energy !
!-------------------------------------!
call wall_time(start_t)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(p,i,a,j,b,n,num,dem1,dem2,reg1,reg2) &
!$OMP SHARED(nC,nO,nOrb,nR,nS,eta,ERI,eQP,eh_trip_rho,eh_trip_Om,SigC,Z)
!$OMP DO COLLAPSE(2)
do p=nC+1,nOrb-nR
! do i=nC+1,nO
! do a=nO+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nOrb-nR
! do n=1,nS
! !3h2p
! do j=nC+1,nO
! num = ( + 1.5d0*ERI(p,a,i,j))* &
! eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n)
do n=1,nS
!3h2p
do j=nC+1,nO
num = ( + 1.5d0*ERI(p,a,i,j))* &
eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n)
! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
! dem2 = eQP(p) - eQP(j) + eh_trip_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
dem2 = eQP(p) - eQP(j) + eh_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( - 1.5d0*ERI(p,a,i,j))* &
! eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n)
num = ( - 1.5d0*ERI(p,a,i,j))* &
eh_trip_rho(a,i,n) * eh_trip_rho(j,p,n)
! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( + 1.5d0*ERI(p,i,a,j)) * &
! eh_trip_rho(j,p,n) * eh_trip_rho(i,a,n)
num = ( + 1.5d0*ERI(p,i,a,j)) * &
eh_trip_rho(j,p,n) * eh_trip_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
! dem2 = eQP(p) - eQP(j) + eh_trip_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
dem2 = eQP(p) - eQP(j) + eh_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( + 1.5d0*ERI(p,a,i,j))* &
! eh_trip_rho(p,j,n) * eh_trip_rho(i,a,n)
num = ( + 1.5d0*ERI(p,a,i,j))* &
eh_trip_rho(p,j,n) * eh_trip_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
! dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
dem2 = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! end do ! j
! !3p2h
! do b=nO+1,nOrb-nR
! num = ( + 1.5d0*ERI(p,a,i,b)) * &
! eh_trip_rho(p,b,n) * eh_trip_rho(i,a,n)
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = ( + 1.5d0*ERI(p,a,i,b)) * &
eh_trip_rho(p,b,n) * eh_trip_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
! dem2 = eQP(p) - eQP(b) - eh_trip_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
dem2 = eQP(p) - eQP(b) - eh_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( + 1.5d0*ERI(p,i,a,b)) * &
! eh_trip_rho(b,p,n) * eh_trip_rho(i,a,n)
num = ( + 1.5d0*ERI(p,i,a,b)) * &
eh_trip_rho(b,p,n) * eh_trip_rho(i,a,n)
! dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) + eh_trip_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( + 1.5d0*ERI(p,i,a,b)) * &
! eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n)
num = ( + 1.5d0*ERI(p,i,a,b)) * &
eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n)
! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
! dem2 = eQP(p) - eQP(b) - eh_trip_Om(n)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
dem2 = eQP(p) - eQP(b) - eh_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! num = ( - 1.5d0*ERI(p,i,a,b)) * &
! eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n)
num = ( - 1.5d0*ERI(p,i,a,b)) * &
eh_trip_rho(p,b,n) * eh_trip_rho(a,i,n)
! dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
! dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
! reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
! reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
dem1 = eQP(a) - eQP(i) - eh_trip_Om(n)
dem2 = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
! end do ! b
end do ! b
! end do ! n
end do ! n
! end do ! a
! end do ! i
end do ! a
end do ! i
! end do ! p
! !$OMP END DO
! !$OMP END PARALLEL
! call wall_time(end_t)
! t = end_t - start_t
end do ! p
!$OMP END DO
!$OMP END PARALLEL
call wall_time(end_t)
t = end_t - start_t
! write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh self-energy =',t,' seconds'
! write(*,*)
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet eh self-energy =',t,' seconds'
write(*,*)
!-------------------------------------!
! singlet pp part of the self-energy !