10
1
mirror of https://github.com/pfloos/quack synced 2025-05-07 07:35:02 +02:00

restricted self energy done

This commit is contained in:
Antoine Marie 2025-04-03 15:50:31 +02:00
parent bf35de22d7
commit 992c65b497
3 changed files with 398 additions and 151 deletions

View File

@ -81,126 +81,126 @@ 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)
! !$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
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
! 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(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
do n=1,nS
!3h2p
do j=nC+1,nO
num = ERI(p,a,j,i) * &
(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(a,p,n) * eh_rho(i,j,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) * &
! !(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
!num = ERI(p,a,j,i) * &
!(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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,a,j,i) * &
! !(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,p,n))
!num = ERI(p,a,j,i) * &
!(eh_rho(j,p,n) * eh_rho(i,a,n) - eh_rho(j,a,n) * eh_rho(i,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) * &
! (eh_rho(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,p,n))
num = ERI(p,i,j,a) * &
(eh_rho(j,p,n) * eh_rho(a,i,n) - eh_rho(i,p,n) * eh_rho(a,j,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)
! end do ! j
! !3p2h
! do b=nO+1,nOrb-nR
! num = ERI(p,a,b,i) * &
! (eh_rho(b,p,n) * eh_rho(i,a,n) - eh_rho(b,a,n) * eh_rho(i,p,n))
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = ERI(p,a,b,i) * &
(eh_rho(b,p,n) * eh_rho(i,a,n) - eh_rho(a,p,n) * eh_rho(i,b,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(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
num = ERI(p,i,b,a) * &
(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(i,p,n) * eh_rho(a,b,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) * &
! !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
!num = ERI(p,i,b,a) * &
!(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,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,i,b,a) * &
! !(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,n))
!num = ERI(p,i,b,a) * &
!(eh_rho(b,p,n) * eh_rho(a,i,n) - eh_rho(b,i,n) * eh_rho(a,p,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)
! 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 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 !
!-----------------------------!

View File

@ -87,6 +87,252 @@ subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds'
write(*,*)
!-------------------------------------!
! 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
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_sing_rho(j,p,n) * eh_sing_rho(i,a,n) - 0.5d0*eh_sing_rho(a,p,n) * eh_sing_rho(i,j,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))
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(j,p,n) * eh_rho(i,a,n) - 0.5d0*eh_rho(j,a,n) * eh_rho(i,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))
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(j,p,n) * eh_rho(i,a,n) - 0.5d0*eh_rho(j,a,n) * eh_rho(i,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))
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_sing_rho(j,p,n) * eh_sing_rho(a,i,n) - 0.5d0*eh_sing_rho(i,p,n) * eh_sing_rho(a,j,n))
!dem1 = eQP(a) - eQP(i) + eh_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)
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = ERI(p,a,b,i) * &
(eh_sing_rho(b,p,n) * eh_sing_rho(i,a,n) - 0.5d0*eh_sing_rho(a,p,n) * eh_sing_rho(i,b,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))
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_sing_rho(b,p,n) * eh_sing_rho(a,i,n) - 0.5d0*eh_sing_rho(i,p,n) * eh_sing_rho(a,b,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)
!num = ERI(p,i,b,a) * &
!(eh_rho(b,p,n) * eh_rho(a,i,n) - 0.5d0*eh_rho(b,i,n) * eh_rho(a,p,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))
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(b,p,n) * eh_rho(a,i,n) - 0.5d0*eh_rho(b,i,n) * eh_rho(a,p,n))
!dem1 = eQP(a) - eQP(i) - eh_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)
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(*,*)
!-------------------------------------!
! 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 n=1,nS
!3h2p
do j=nC+1,nO
num = ERI(p,a,j,i) * &
(- 1.5d0*eh_trip_rho(a,p,n) * eh_trip_rho(i,j,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))
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) * &
!(- 1.5d0*eh_rho(j,a,n) * eh_rho(i,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))
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) * &
!(- 1.5d0*eh_rho(j,a,n) * eh_rho(i,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))
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) * &
(- 1.5d0*eh_trip_rho(i,p,n) * eh_trip_rho(a,j,n))
!dem1 = eQP(a) - eQP(i) + eh_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)
end do ! j
!3p2h
do b=nO+1,nOrb-nR
num = ERI(p,a,b,i) * &
(- 1.5d0*eh_trip_rho(a,p,n) * eh_trip_rho(i,b,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))
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) * &
(- 1.5d0*eh_trip_rho(i,p,n) * eh_trip_rho(a,b,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)
!num = ERI(p,i,b,a) * &
!(eh_rho(b,p,n) * eh_rho(a,i,n) - 0.5d0*eh_rho(b,i,n) * eh_rho(a,p,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))
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(b,p,n) * eh_rho(a,i,n) - 0.5d0*eh_rho(b,i,n) * eh_rho(a,p,n))
!dem1 = eQP(a) - eQP(i) - eh_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)
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(*,*)
!-------------------------------------!
! singlet pp part of the self-energy !
!-------------------------------------!
call wall_time(start_t)

View File

@ -38,12 +38,12 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) &
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) - 0d0*1.5d0*eh_trip_Phi(q,j,b,p) &
+ 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*1.5d0*pp_trip_Phi(q,j,p,b)) * X &
+ (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) &
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) - 0d0*1.5d0*eh_trip_Phi(q,b,j,p) &
+ 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*1.5d0*pp_trip_Phi(q,b,p,j)) * Y
rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(q,j,p,b) - ERI(q,j,b,p))*X ! &
!(- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) - 0d0*1.5d0*eh_trip_Phi(q,j,b,p) &
!+ 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*1.5d0*pp_trip_Phi(q,j,p,b)) * X &
!+ (2d0*ERI(q,b,p,j) - ERI(q,b,j,p))*Y &
!(- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) - 0d0*1.5d0*eh_trip_Phi(q,b,j,p) &
!+ 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*1.5d0*pp_trip_Phi(q,b,p,j)) * Y
end do
@ -97,12 +97,13 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr
X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb))
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
rho(p,q,ia) = rho(p,q,ia) + (- ERI(q,j,b,p) &
- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) + 0d0*0.5d0*eh_trip_Phi(q,j,b,p) &
- 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*0.5d0*pp_trip_Phi(q,j,p,b)) * X &
+ (- ERI(q,b,j,p) &
- 0d0*0.5d0*eh_sing_Phi(q,b,j,p) + 0d0*0.5d0*eh_trip_Phi(q,b,j,p) &
- 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*0.5d0*pp_trip_Phi(q,b,p,j)) * Y
rho(p,q,ia) = rho(p,q,ia) - ERI(q,j,b,p) * X ! &
! (- 0d0*0.5d0*eh_sing_Phi(q,j,b,p) + 0d0*0.5d0*eh_trip_Phi(q,j,b,p) &
! - 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*0.5d0*pp_trip_Phi(q,j,p,b)) * X &
! + (- ERI(q,b,j,p) &
! - 0d0*0.5d0*eh_sing_Phi(q,b,j,p) + 0d0*0.5d0*eh_trip_Phi(q,b,j,p) &
! - 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*0.5d0*pp_trip_Phi(q,b,p,j)) * Y
end do
end do
@ -173,26 +174,26 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do d = c, nOrb-nR
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) + ERI(p,q,d,c)) &
!+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
!+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
*X1(cd,ab)/sqrt(1d0 + Kronecker_delta(c,d))
end do ! d
end do ! c
kl = 0
do k = nC+1, nO
do l = k, nO
! kl = 0
! do k = nC+1, nO
! do l = k, nO
kl = kl + 1
! kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
*Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l))
end do ! l
end do ! k
! rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
! *Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l))
! end do ! l
! end do ! k
end do ! b
end do ! a
@ -201,26 +202,26 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do j = i, nO
ij = ij + 1
cd = 0
do c = nO+1, nOrb-nR
do d = c, nOrb-nR
cd = cd + 1
! cd = 0
! do c = nO+1, nOrb-nR
! do d = c, nOrb-nR
! cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
*X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d))
end do ! d
end do ! c
! rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) + ERI(p,q,d,c) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))&
! *X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d))
! end do ! d
! end do ! c
kl = 0
do k = nC+1, nO
do l = k, nO
kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) + ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) + ERI(p,q,l,k)) &
!+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) &
!+ 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))&
*Y2(kl,ij)/sqrt(1d0 + Kronecker_delta(k,l))
end do ! l
end do ! k
@ -290,24 +291,24 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do d = c+1, nOrb-nR
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab)
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)! &
!+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
!- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab)
end do ! d
end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
! kl = 0
! do k = nC+1, nO
! do l = k+1, nO
kl = kl + 1
! kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab)
end do ! l
end do ! k
! rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
! - 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab)
! end do ! l
! end do ! k
end do ! b
end do ! a
@ -316,25 +317,25 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi,
do j = i+1, nO
ij = ij + 1
cd = 0
do c = nO+1, nOrb-nR
do d = c+1, nOrb-nR
cd = cd + 1
! cd = 0
! do c = nO+1, nOrb-nR
! do d = c+1, nOrb-nR
! cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
- 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij)
end do ! d
end do ! c
! rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c) &
! + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) &
! - 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij)
! end do ! d
! end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k) &
+ 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
- 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij)
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)! &
! + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) &
! - 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij)
end do ! l
end do ! k