10
1
mirror of https://github.com/pfloos/quack synced 2025-05-06 15:14:55 +02:00

GW alternatif

This commit is contained in:
Antoine Marie 2025-04-23 14:19:29 +02:00
parent 921902962b
commit 09051226de
5 changed files with 289 additions and 283 deletions

View File

@ -43,129 +43,129 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh
! Occupied part of the Tpp self-energy !
!--------------------------------------!
do p=nC+1,nBas-nR
do i=nC+1,nO
! do p=nC+1,nBas-nR
! do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Om1(cd)
num = rho1(p,i,cd)**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 cd=1,nVV
! eps = e(p) + e(i) - Om1(cd)
! num = rho1(p,i,cd)**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
end do
end do
! end do
! end do
!------------------------------------------!
! Virtual part of the T-matrix self-energy !
!------------------------------------------!
! !------------------------------------------!
! ! Virtual part of the T-matrix self-energy !
! !------------------------------------------!
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
! do p=nC+1,nBas-nR
! do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Om2(kl)
num = rho2(p,a,kl)**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 kl=1,nOO
! eps = e(p) + e(a) - Om2(kl)
! num = rho2(p,a,kl)**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
end do
end do
! 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
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
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
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
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,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
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
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
! 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
! 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
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
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
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,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
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
! 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
! end do
end do
end do
end do
!-------------------------------------!
! Galitskii-Migdal correlation energy !

View File

@ -119,7 +119,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
if(doSRG) then
call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
else
call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z,ERI)
end if
!-----------------------------------!

View File

@ -1,4 +1,4 @@
subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z,ERI)
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
@ -17,11 +17,12 @@ subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,a,p,m
double precision :: num,eps
integer :: i,j,a,b,p,m
double precision :: num,eps,dem1,dem2
! Output variables
@ -38,36 +39,118 @@ subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
! ! Occupied part of the correlation self-energy
! do p=nC+1,nBas-nR
! do i=nC+1,nO
! do m=1,nS
! eps = e(p) - e(i) + Om(m)
! num = rho(p,i,m)**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
! end do
! end do
! ! Virtual part of the correlation self-energy
! do p=nC+1,nBas-nR
! do a=nO+1,nBas-nR
! do m=1,nS
! eps = e(p) - e(a) - Om(m)
! num = rho(p,a,m)**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
! end do
! end do
!-----------------------------------------------!
! Testing another way to compute GT self-energy !
!-----------------------------------------------!
do p=nC+1,nBas-nR
do i=nC+1,nO
do m=1,nS
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = e(p) - e(i) + Om(m)
num = rho(p,i,m)**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
eps = e(p) + e(a) - e(i) - e(j)
num = ERI(p,a,i,j)**2
end do
end do
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,nS
num = - ERI(p,i,j,a) * rho(i,a,m) * rho(j,p,m)
dem1 = e(p) + e(a) - e(i) - e(j)
dem2 = e(p) - e(j) + Om(m)
Sig(p) = Sig(p) + num/dem1/dem2
Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2
num = - ERI(p,i,j,a) * rho(i,a,m) * rho(j,p,m)
dem1 = e(p) + e(a) - e(i) - e(j)
dem2 = e(a) - e(i) + Om(m)
Sig(p) = Sig(p) + num/dem1/dem2
Z(p) = Z(p) - num/dem1/dem1/dem2
num = - ERI(p,a,j,i) * rho(a,i,m) * rho(j,p,m)
dem1 = e(p) - e(j) + Om(m)
dem2 = e(a) - e(i) + Om(m)
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
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do m=1,nS
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do i=nC+1,nO
eps = e(p) - e(a) - Om(m)
num = rho(p,a,m)**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
eps = e(p) + e(i) - e(a) - e(b)
num = ERI(p,i,a,b)**2
end do
end do
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,nS
num = ERI(p,a,b,i) * rho(a,i,m) * rho(b,p,m)
dem1 = e(p) + e(i) - e(a) - e(b)
dem2 = e(p) - e(b) - Om(m)
Sig(p) = Sig(p) + num/dem1/dem2
Z(p) = Z(p) - num/dem1/dem1/dem2 - num/dem1/dem2/dem2
num = - ERI(p,a,b,i) * rho(a,i,m) * rho(b,p,m)
dem1 = e(p) + e(i) - e(a) - e(b)
dem2 = e(a) - e(i) + Om(m)
Sig(p) = Sig(p) + num/dem1/dem2
Z(p) = Z(p) - num/dem1/dem1/dem2
num = - ERI(p,i,b,a) * rho(i,a,m) * rho(b,p,m)
dem1 = e(p) - e(b) - Om(m)
dem2 = e(a) - e(i) + Om(m)
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
EcGM = 0d0

View File

@ -55,6 +55,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
eps = eQP(p) + eQP(a) - eQP(i) - eQP(j)
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
num = 0.5d0*(ERI(p,a,j,i) - ERI(p,a,i,j))**2
! num = ERI(p,a,j,i)**2
SigC(p) = SigC(p) + num*reg/eps
Z(p) = Z(p) - num*reg/eps**2
@ -70,6 +71,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
eps = eQP(p) + eQP(i) - eQP(a) - eQP(b)
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
num = 0.5d0*(ERI(p,i,b,a) - ERI(p,i,a,b))**2
! num = ERI(p,i,b,a)**2
SigC(p) = SigC(p) + num*reg/eps
Z(p) = Z(p) - num*reg/eps**2
@ -89,151 +91,108 @@ 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
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(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))
num = - (ERI(p,a,j,i) - ERI(p,a,i,j)) * 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)
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 = - 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))
num = - (ERI(p,i,j,a) - ERI(p,i,a,j)) * 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)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/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))
! ! 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)
SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
Z(p) = Z(p) - num * (reg1/dem1) * (reg2/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))
num = - (ERI(p,a,j,i) - ERI(p,a,i,j)) * 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)
! ! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! ! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/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))
! 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)
! 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) - ERI(p,i,a,b)) * 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(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))
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) &
! - num * (reg1/dem1/dem1) * (reg2/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,a,b,i) - ERI(p,a,i,b)) * 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(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) * 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) - ERI(p,i,a,b)) * 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))
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,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)
end do ! b
! 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))
end do ! n
! SigC(p) = SigC(p) + num * (reg1/dem1) * (reg2/dem2)
! Z(p) = Z(p) - num * (reg1/dem1) * (reg2/dem2/dem2)
end do ! a
end do ! i
! 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)
end do ! p
!$OMP END DO
!$OMP END PARALLEL
! 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))
call wall_time(end_t)
t = end_t - start_t
! 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 ! a
! end do ! i
! 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(*,*)
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds'
write(*,*)
!-----------------------------!
! pp part of the self-energy !
@ -290,24 +249,6 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
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
end do ! j
@ -365,24 +306,6 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
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))
! 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
end do ! b

View File

@ -36,15 +36,15 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
rho(p,q,ia) = rho(p,q,ia) &
+ (ERI(q,j,p,b) - ERI(q,j,b,p)) * X !&
+ (ERI(q,j,p,b) - ERI(q,j,b,p)) * X &
!- (eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X &
!+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y &
+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y !&
!- (eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y
rho(p,q,nS+ia) = rho(p,q,nS+ia) &
+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * X !&
+ (ERI(q,j,p,b) - ERI(q,b,j,p)) * X &
!- (eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X &
!+ (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y &
+ (ERI(q,b,p,j) - ERI(q,j,b,p)) * Y !&
!- (eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y
end do