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

ok with pp self-energy, eh self-energy still giving crazy results

Browse Source
This commit is contained in:
Antoine Marie 2025-04-15 15:55:24 +02:00
parent 095695ab6b
commit 1d316f0fd2
7 changed files with 348 additions and 346 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Parquet/GParquet.f90
View File

@ -8,8 +8,8 @@ subroutine GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
! Hard-coded parameters
logical :: print_phLR = .false.
logical :: print_ppLR = .false.
logical :: print_phLR = .true.
logical :: print_ppLR = .true.
! Input variables

212
src/Parquet/G_Parquet_self_energy.f90
View File

@ -89,129 +89,129 @@ 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) + ERI(p,a,i,j))* &
eh_rho(a,i,n) * eh_rho(j,p,n)
! do n=1,nS
! !3h2p
! do j=nC+1,nO
! num = ( - ERI(p,a,j,i) + 0d0*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(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)
! 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) - ERI(p,a,i,j))* &
eh_rho(a,i,n) * eh_rho(j,p,n)
! num = (ERI(p,a,j,i) - 0d0*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(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)
num = (- ERI(p,i,j,a) + ERI(p,i,a,j)) * &
eh_rho(j,p,n) * eh_rho(i,a,n)
! num = (- ERI(p,i,j,a) + 0d0*ERI(p,i,a,j)) * &
! eh_rho(j,p,n) * eh_rho(i,a,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(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)
! 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) + ERI(p,a,i,j))* &
eh_rho(p,j,n) * eh_rho(i,a,n)
! num = (- ERI(p,a,j,i) + 0d0*ERI(p,a,i,j))* &
! eh_rho(p,j,n) * eh_rho(i,a,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,a,b,i) + ERI(p,a,i,b)) * &
eh_rho(p,b,n) * eh_rho(i,a,n)
! end do ! j
! !3p2h
! do b=nO+1,nOrb-nR
! num = (- ERI(p,a,b,i) + 0d0*ERI(p,a,i,b)) * &
! eh_rho(p,b,n) * eh_rho(i,a,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(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,i,b,a) + ERI(p,i,a,b)) * &
eh_rho(b,p,n) * eh_rho(i,a,n)
! num = (- ERI(p,i,b,a) + 0d0*ERI(p,i,a,b)) * &
! eh_rho(b,p,n) * eh_rho(i,a,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(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) + ERI(p,i,a,b)) * &
eh_rho(p,b,n) * eh_rho(a,i,n)
! num = (- ERI(p,i,b,a) + 0d0*ERI(p,i,a,b)) * &
! eh_rho(p,b,n) * eh_rho(a,i,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(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,i,b,a) - ERI(p,i,a,b)) * &
eh_rho(p,b,n) * eh_rho(a,i,n)
! num = (ERI(p,i,b,a) - 0d0*ERI(p,i,a,b)) * &
! eh_rho(p,b,n) * eh_rho(a,i,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(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
! 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 building eh self-energy =',t,' seconds'
write(*,*)
! 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(*,*)
!-----------------------------!
! pp part of the self-energy !
!-----------------------------!
@ -229,7 +229,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
do n=1,nVV
! 4h1p
do k=nC+1,nO
num = 2d0 * ERI(p,k,i,j) * ee_rho(i,j,n) * ee_rho(p,k,n)
num = - ERI(p,k,i,j) * ee_rho(i,j,n) * ee_rho(p,k,n)
dem1 = ee_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(k) - ee_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -242,7 +242,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
! 3h2p
do c=nO+1,nOrb-nR
num = 2d0 * ERI(p,c,i,j) * ee_rho(i,j,n) * ee_rho(p,c,n)
num = - ERI(p,c,i,j) * ee_rho(i,j,n) * ee_rho(p,c,n)
dem1 = ee_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -257,7 +257,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
! 3h2p
do c=nO+1,nOrb-nR
num = 2d0 * ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
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))
@ -266,14 +266,14 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = 2d0 * ERI(p,c,i,j) * hh_rho(i,j,n) * hh_rho(p,c,n)
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)
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
@ -294,27 +294,27 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
! 4p1h
do c=nO+1,nOrb-nR
num = 2d0 * ERI(p,c,a,b) * hh_rho(a,b,n) * hh_rho(p,c,n)
num = ERI(p,c,a,b) * hh_rho(a,b,n) * hh_rho(p,c,n)
dem1 = hh_Om(n) - eQP(a) - eQP(b)
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)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! c
! 3p2h
do k=nC+1,nO
num = 2d0 * ERI(p,k,a,b) * hh_rho(a,b,n) * hh_rho(p,k,n)
num = ERI(p,k,a,b) * hh_rho(a,b,n) * hh_rho(p,k,n)
dem1 = hh_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)
end do ! k
end do ! n
@ -322,16 +322,16 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
! 3p2h
do k=nC+1,nO
num = 2d0 * ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
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))
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 = 2d0 * ERI(p,k,a,b) * ee_rho(a,b,n) * ee_rho(p,k,n)
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))

4
src/Parquet/G_pp_Phi.f90
View File

@ -31,12 +31,12 @@ subroutine G_pp_Phi(nOrb,nC,nR,nOO,nVV,ee_Om,ee_rho,hh_Om,hh_rho,pp_Phi)
do n=1,nVV
pp_Phi(p,q,r,s) = pp_Phi(p,q,r,s) &
+ 2d0 * ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n)
- ee_rho(p,q,n)*ee_rho(r,s,n)/ee_Om(n)
end do
do n=1,nOO
pp_Phi(p,q,r,s) = pp_Phi(p,q,r,s) &
- 2d0 * hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n)
+ hh_rho(p,q,n)*hh_rho(r,s,n)/hh_Om(n)
end do
enddo

4
src/Parquet/RParquet.f90
View File

@ -8,8 +8,8 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
! Hard-coded parameters
logical :: print_phLR = .false.
logical :: print_ppLR = .false.
logical :: print_phLR = .true.
logical :: print_ppLR = .true.
! Input variables

462
src/Parquet/R_Parquet_self_energy.f90
View File

@ -89,249 +89,250 @@ 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 !
!-------------------------------------!
@ -347,7 +348,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
do n=1,nVVs
! 4h1p
do k=nC+1,nO
num = ERI(p,k,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,k,n)
num = - 0.5d0 * ERI(p,k,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,k,n)
dem1 = ee_sing_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(k) - ee_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -360,10 +361,10 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3h2p
do c=nO+1,nOrb-nR
num = ERI(p,c,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,c,n)
!dem1 = ee_Om(n) - eQP(i) - eQP(j)
num = - 0.5d0*ERI(p,c,i,j) * ee_sing_rho(i,j,n) * ee_sing_rho(p,c,n)
dem1 = ee_sing_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
@ -375,7 +376,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3h2p
do c=nO+1,nOrb-nR
num = ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n)
num = - 0.5d0*ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n)
dem1 = hh_sing_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - hh_sing_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -384,14 +385,14 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
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_sing_rho(i,j,n) * hh_sing_rho(p,c,n)
!dem1 = hh_Om(n) - eQP(i) - eQP(j)
num = 0.5d0*ERI(p,c,i,j) * hh_sing_rho(i,j,n) * hh_sing_rho(p,c,n)
dem1 = hh_sing_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
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 ! c
end do ! n
@ -412,27 +413,27 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 4p1h
do c=nO+1,nOrb-nR
num = ERI(p,c,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,c,n)
num = 0.5d0*ERI(p,c,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,c,n)
dem1 = hh_sing_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(c) - hh_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)
end do ! c
! 3p2h
do k=nC+1,nO
num = ERI(p,k,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,k,n)
!dem1 = hh_Om(n) - eQP(a) - eQP(b)
num = 0.5d0*ERI(p,k,a,b) * hh_sing_rho(a,b,n) * hh_sing_rho(p,k,n)
dem1 = hh_sing_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
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 ! k
end do ! n
@ -440,19 +441,19 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3p2h
do k=nC+1,nO
num = ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n)
num = 0.5d0*ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n)
dem1 = ee_sing_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - ee_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,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n)
!dem1 = ee_Om(n) - eQP(a) - eQP(b)
num = - 0.5d0*ERI(p,k,a,b) * ee_sing_rho(a,b,n) * ee_sing_rho(p,k,n)
dem1 = ee_sing_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
@ -486,7 +487,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
do n=1,nVVt
! 4h1p
do k=nC+1,nO
num = 3d0*ERI(p,k,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,k,n)
num = - 1.5d0 * ERI(p,k,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,k,n)
dem1 = ee_trip_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(k) - ee_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -499,10 +500,10 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3h2p
do c=nO+1,nOrb-nR
num = 3d0*ERI(p,c,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,c,n)
!dem1 = ee_Om(n) - eQP(i) - eQP(j)
num = - 1.5d0 * ERI(p,c,i,j) * ee_trip_rho(i,j,n) * ee_trip_rho(p,c,n)
dem1 = ee_trip_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
@ -514,7 +515,7 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3h2p
do c=nO+1,nOrb-nR
num = 3d0*ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n)
num = - 1.5d0 * ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n)
dem1 = hh_trip_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - hh_trip_Om(n)
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
@ -523,14 +524,14 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
num = 3d0*ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n)
!dem1 = hh_Om(n) - eQP(i) - eQP(j)
num = 1.5d0 * ERI(p,c,i,j) * hh_trip_rho(i,j,n) * hh_trip_rho(p,c,n)
dem1 = hh_trip_Om(n) - eQP(i) - eQP(j)
dem2 = eQP(p) + eQP(c) - eQP(i) - eQP(j)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
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 ! c
end do ! n
@ -551,27 +552,27 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 4p1h
do c=nO+1,nOrb-nR
num = 3d0*ERI(p,c,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,c,n)
num = 1.5d0 * ERI(p,c,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,c,n)
dem1 = hh_trip_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(c) - hh_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)
end do ! c
! 3p2h
do k=nC+1,nO
num = 3d0*ERI(p,k,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,k,n)
!dem1 = hh_Om(n) - eQP(a) - eQP(b)
num = 1.5d0 * ERI(p,k,a,b) * hh_trip_rho(a,b,n) * hh_trip_rho(p,k,n)
dem1 = hh_trip_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
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 ! k
end do ! n
@ -579,19 +580,19 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
! 3p2h
do k=nC+1,nO
num = 3d0*ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n)
num = 1.5d0 * ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n)
dem1 = ee_trip_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - ee_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 = 3d0*ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n)
!dem1 = ee_Om(n) - eQP(a) - eQP(b)
num = - 1.5d0 * ERI(p,k,a,b) * ee_trip_rho(a,b,n) * ee_trip_rho(p,k,n)
dem1 = ee_trip_Om(n) - eQP(a) - eQP(b)
dem2 = eQP(p) + eQP(k) - eQP(a) - eQP(b)
!reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg1 = (1d0 - exp(- 2d0 * eta * dem1 * dem1))
reg2 = (1d0 - exp(- 2d0 * eta * dem2 * dem2))
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
@ -608,7 +609,8 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
call wall_time(end_t)
t = end_t - start_t
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for triplet pp self-energy =',t,' seconds'
write(*,*)
write(*,*)
!-----------------------------!
! Renormalization factor !
!-----------------------------!

4
src/Parquet/R_pp_singlet_Phi.f90
View File

@ -32,12 +32,12 @@ subroutine R_pp_singlet_Phi(nOrb,nC,nR,nOO,nVV,ee_sing_Om,ee_sing_rho,hh_sing_Om
do n=1,nVV
pp_sing_Phi(p,q,r,s) = pp_sing_Phi(p,q,r,s) &
+ 2d0 * ee_sing_rho(p,q,n)*ee_sing_rho(r,s,n)/ee_sing_Om(n)
- ee_sing_rho(p,q,n)*ee_sing_rho(r,s,n)/ee_sing_Om(n)
end do
do n=1,nOO
pp_sing_Phi(p,q,r,s) = pp_sing_Phi(p,q,r,s) &
- 2d0 * hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n)
+ hh_sing_rho(p,q,n)*hh_sing_rho(r,s,n)/hh_sing_Om(n)
end do
enddo

4
src/Parquet/R_pp_triplet_Phi.f90
View File

@ -32,12 +32,12 @@ subroutine R_pp_triplet_Phi(nOrb,nC,nR,nOO,nVV,ee_trip_Om,ee_trip_rho,hh_trip_Om
do n=1,nVV
pp_trip_Phi(p,q,r,s) = pp_trip_Phi(p,q,r,s) &
+ 2d0 * ee_trip_rho(p,q,n)*ee_trip_rho(r,s,n)/ee_trip_Om(n)
- ee_trip_rho(p,q,n)*ee_trip_rho(r,s,n)/ee_trip_Om(n)
end do
do n=1,nOO
pp_trip_Phi(p,q,r,s) = pp_trip_Phi(p,q,r,s) &
- 2d0 * hh_trip_rho(p,q,n)*hh_trip_rho(r,s,n)/hh_trip_Om(n)
+ hh_trip_rho(p,q,n)*hh_trip_rho(r,s,n)/hh_trip_Om(n)
end do
enddo