mirror of
https://github.com/pfloos/quack
synced 2025-05-06 23:24:58 +02:00
Merge branch 'master' of github.com:pfloos/QuAcK
This commit is contained in:
commit
2053384ff3
@ -554,8 +554,8 @@ subroutine GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
|
|||||||
eOld(:) = eQP(:)
|
eOld(:) = eQP(:)
|
||||||
|
|
||||||
! Print for one-body part
|
! Print for one-body part
|
||||||
|
|
||||||
call print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp)
|
call G_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp)
|
||||||
|
|
||||||
deallocate(eQPlin,Z,SigC)
|
deallocate(eQPlin,Z,SigC)
|
||||||
|
|
||||||
|
|||||||
356
src/Parquet/G_Parquet_self_energy.f90
Normal file
356
src/Parquet/G_Parquet_self_energy.f90
Normal file
@ -0,0 +1,356 @@
|
|||||||
|
subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
|
||||||
|
eh_rho,eh_Om,ee_rho,ee_Om,hh_rho,hh_Om,EcGM,SigC,Z)
|
||||||
|
|
||||||
|
! Compute correlation part of the self-energy coming from irreducible vertices contribution
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
integer,intent(in) :: nOrb
|
||||||
|
integer,intent(in) :: nC, nO, nV, nR
|
||||||
|
integer,intent(in) :: nS, nOO, nVV
|
||||||
|
double precision,intent(in) :: eQP(nOrb)
|
||||||
|
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||||
|
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
|
||||||
|
double precision,intent(in) :: eh_Om(nS)
|
||||||
|
double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV)
|
||||||
|
double precision,intent(in) :: ee_Om(nVV)
|
||||||
|
double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO)
|
||||||
|
double precision,intent(in) :: hh_Om(nOO)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
integer :: i,j,k,a,b,c
|
||||||
|
integer :: p,n
|
||||||
|
double precision :: eps,dem1,dem2,reg,reg1,reg2
|
||||||
|
double precision :: num
|
||||||
|
double precision :: start_t,end_t,t
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
double precision,intent(out) :: SigC(nOrb)
|
||||||
|
double precision,intent(out) :: Z(nOrb)
|
||||||
|
double precision,intent(out) :: EcGM
|
||||||
|
|
||||||
|
! Initialize
|
||||||
|
|
||||||
|
SigC(:) = 0d0
|
||||||
|
Z(:) = 0d0
|
||||||
|
EcGM = 0d0
|
||||||
|
|
||||||
|
!-----------------------------!
|
||||||
|
! GF2 part of the self-energy !
|
||||||
|
!-----------------------------!
|
||||||
|
call wall_time(start_t)
|
||||||
|
do p=nC+1,nOrb-nR
|
||||||
|
! 2h1p sum
|
||||||
|
do i=nC+1,nO
|
||||||
|
do j=nC+1,nO
|
||||||
|
do a=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
|
SigC(p) = SigC(p) + num*reg/eps
|
||||||
|
Z(p) = Z(p) - num*reg/eps**2
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 2p1h sum
|
||||||
|
do i=nC+1,nO
|
||||||
|
do a=nO+1,nOrb-nR
|
||||||
|
do b=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
|
SigC(p) = SigC(p) + num*reg/eps
|
||||||
|
Z(p) = Z(p) - num*reg/eps**2
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
call wall_time(end_t)
|
||||||
|
t = end_t - start_t
|
||||||
|
|
||||||
|
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds'
|
||||||
|
write(*,*)
|
||||||
|
! !-----------------------------!
|
||||||
|
! ! 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
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
|
||||||
|
! 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(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))
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
|
||||||
|
! 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,i,j,a) * &
|
||||||
|
! (eh_rho(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,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)
|
||||||
|
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
|
||||||
|
! 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)
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
|
||||||
|
! 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) - 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))
|
||||||
|
|
||||||
|
! 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(*,*)
|
||||||
|
! !-----------------------------!
|
||||||
|
! ! pp part of the self-energy !
|
||||||
|
! !-----------------------------!
|
||||||
|
! call wall_time(start_t)
|
||||||
|
! !$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
! !$OMP PRIVATE(p,i,j,k,c,n,num,dem1,dem2,reg1,reg2) &
|
||||||
|
! !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z)
|
||||||
|
! !$OMP DO COLLAPSE(2)
|
||||||
|
! do p=nC+1,nOrb-nR
|
||||||
|
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do n=1,nVV
|
||||||
|
! ! 4h1p
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! 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))
|
||||||
|
! 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 ! k
|
||||||
|
! ! 3h2p
|
||||||
|
! do c=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
! 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))
|
||||||
|
! 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 ! a
|
||||||
|
! end do ! n
|
||||||
|
! do n=1,nOO
|
||||||
|
! ! 3h2p
|
||||||
|
! 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)
|
||||||
|
! 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
|
||||||
|
! end do ! i
|
||||||
|
|
||||||
|
! end do ! p
|
||||||
|
! !$OMP END DO
|
||||||
|
! !$OMP END PARALLEL
|
||||||
|
! !$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
! !$OMP PRIVATE(p,k,a,b,c,n,num,dem1,dem2,reg1,reg2) &
|
||||||
|
! !$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z)
|
||||||
|
! !$OMP DO COLLAPSE(2)
|
||||||
|
! do p=nC+1,nOrb-nR
|
||||||
|
! do a=nO+1,nOrb-nR
|
||||||
|
! do b=nO+1,nOrb-nR
|
||||||
|
! do n=1,nOO
|
||||||
|
! ! 4p1h
|
||||||
|
! do c=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
! 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)
|
||||||
|
|
||||||
|
! end do ! c
|
||||||
|
! ! 3p2h
|
||||||
|
! do k=nC+1,nO
|
||||||
|
|
||||||
|
! 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)
|
||||||
|
|
||||||
|
! end do ! k
|
||||||
|
! end do ! n
|
||||||
|
! do n=1,nVV
|
||||||
|
! ! 3p2h
|
||||||
|
! 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)
|
||||||
|
! 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
|
||||||
|
! end do ! a
|
||||||
|
|
||||||
|
! 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 pp self-energy =',t,' seconds'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
!-----------------------------!
|
||||||
|
! Renormalization factor !
|
||||||
|
!-----------------------------!
|
||||||
|
|
||||||
|
Z(:) = 1d0/(1d0 - Z(:))
|
||||||
|
|
||||||
|
!-------------------------------------!
|
||||||
|
! Galitskii-Migdal correlation energy !
|
||||||
|
!-------------------------------------!
|
||||||
|
|
||||||
|
EcGM = 0d0
|
||||||
|
|
||||||
|
end subroutine
|
||||||
@ -1,349 +0,0 @@
|
|||||||
subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,&
|
|
||||||
eh_rho,eh_Om,ee_rho,ee_Om,hh_rho,hh_Om,EcGM,SigC,Z)
|
|
||||||
|
|
||||||
! Compute correlation part of the self-energy coming from irreducible vertices contribution
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
double precision,intent(in) :: eta
|
|
||||||
integer,intent(in) :: nOrb
|
|
||||||
integer,intent(in) :: nC, nO, nV, nR
|
|
||||||
integer,intent(in) :: nS, nOO, nVV
|
|
||||||
double precision,intent(in) :: eQP(nOrb)
|
|
||||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
|
||||||
double precision,intent(in) :: eh_rho(nOrb,nOrb,nS)
|
|
||||||
double precision,intent(in) :: eh_Om(nS)
|
|
||||||
double precision,intent(in) :: ee_rho(nOrb,nOrb,nVV)
|
|
||||||
double precision,intent(in) :: ee_Om(nVV)
|
|
||||||
double precision,intent(in) :: hh_rho(nOrb,nOrb,nOO)
|
|
||||||
double precision,intent(in) :: hh_Om(nOO)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
integer :: i,j,k,a,b,c
|
|
||||||
integer :: p,n
|
|
||||||
double precision :: eps,dem1,dem2,reg,reg1,reg2
|
|
||||||
double precision :: num
|
|
||||||
double precision :: start_t,end_t,t
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
double precision,intent(out) :: SigC(nOrb)
|
|
||||||
double precision,intent(out) :: Z(nOrb)
|
|
||||||
double precision,intent(out) :: EcGM
|
|
||||||
|
|
||||||
! Initialize
|
|
||||||
|
|
||||||
SigC(:) = 0d0
|
|
||||||
Z(:) = 0d0
|
|
||||||
EcGM = 0d0
|
|
||||||
|
|
||||||
!-----------------------------!
|
|
||||||
! GF2 part of the self-energy !
|
|
||||||
!-----------------------------!
|
|
||||||
call wall_time(start_t)
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
! 2h1p sum
|
|
||||||
do i=nC+1,nO
|
|
||||||
do j=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
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
|
|
||||||
|
|
||||||
SigC(p) = SigC(p) + num*reg/eps
|
|
||||||
Z(p) = Z(p) - num*reg/eps**2
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2p1h sum
|
|
||||||
do i=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
do b=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
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
|
|
||||||
|
|
||||||
SigC(p) = SigC(p) + num*reg/eps
|
|
||||||
Z(p) = Z(p) - num*reg/eps**2
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
call wall_time(end_t)
|
|
||||||
t = end_t - start_t
|
|
||||||
|
|
||||||
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds'
|
|
||||||
write(*,*)
|
|
||||||
!-----------------------------!
|
|
||||||
! 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
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
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(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))
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
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,i,j,a) * &
|
|
||||||
(eh_rho(j,p,n) * eh_rho(a,i,n) - eh_rho(j,i,n) * eh_rho(a,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)
|
|
||||||
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
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))
|
|
||||||
|
|
||||||
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) - 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))
|
|
||||||
|
|
||||||
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(*,*)
|
|
||||||
!-----------------------------!
|
|
||||||
! pp part of the self-energy !
|
|
||||||
!-----------------------------!
|
|
||||||
call wall_time(start_t)
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(p,i,j,k,c,n,num,dem1,dem2,reg1,reg2) &
|
|
||||||
!$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z)
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
|
|
||||||
do i=nC+1,nO
|
|
||||||
do j=nC+1,nO
|
|
||||||
do n=1,nVV
|
|
||||||
! 4h1p
|
|
||||||
do k=nC+1,nO
|
|
||||||
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))
|
|
||||||
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 ! k
|
|
||||||
! 3h2p
|
|
||||||
do c=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
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))
|
|
||||||
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 ! a
|
|
||||||
end do ! n
|
|
||||||
do n=1,nOO
|
|
||||||
! 3h2p
|
|
||||||
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)
|
|
||||||
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
|
|
||||||
end do ! i
|
|
||||||
|
|
||||||
end do ! p
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(p,k,a,b,c,n,num,dem1,dem2,reg1,reg2) &
|
|
||||||
!$OMP SHARED(nC,nO,nOrb,nR,nOO,nVV,eta,ERI,eQP,ee_rho,ee_Om,hh_rho,hh_Om,SigC,Z)
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
do b=nO+1,nOrb-nR
|
|
||||||
do n=1,nOO
|
|
||||||
! 4p1h
|
|
||||||
do c=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
end do ! c
|
|
||||||
! 3p2h
|
|
||||||
do k=nC+1,nO
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
end do ! k
|
|
||||||
end do ! n
|
|
||||||
do n=1,nVV
|
|
||||||
! 3p2h
|
|
||||||
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)
|
|
||||||
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
|
|
||||||
end do ! a
|
|
||||||
|
|
||||||
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 pp self-energy =',t,' seconds'
|
|
||||||
write(*,*)
|
|
||||||
!-----------------------------!
|
|
||||||
! Renormalization factor !
|
|
||||||
!-----------------------------!
|
|
||||||
|
|
||||||
Z(:) = 1d0/(1d0 - Z(:))
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp)
|
subroutine G_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp)
|
||||||
|
|
||||||
! Print one-electron energies and other stuff for G0F2
|
! Print one-electron energies and other stuff for G0F2
|
||||||
|
|
||||||
@ -57,7 +57,7 @@ subroutine print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM
|
|||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,'(2X,A60,F15.6,A3)') ' eh-RPA correlation energy = ',Ec_eh,' au'
|
write(*,'(2X,A60,F15.6,A3)') ' eh-RPA correlation energy = ',Ec_eh,' au'
|
||||||
write(*,'(2X,A60,F15.6,A3)') ' pp-RPA correlation energy = ',Ec_pp,' au'
|
write(*,'(2X,A60,F15.6,A3)') ' pp-RPA correlation energy = ',Ec_pp,' au'
|
||||||
write(*,'(2X,A60,F15.6,A3)') '(eh+pp)-RPA correlation energy = ',Ec_pp,' au'
|
!write(*,'(2X,A60,F15.6,A3)') '(eh+pp)-RPA correlation energy = ',Ec_pp,' au'
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
@ -600,7 +600,7 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
|
|||||||
write(*,*)' Two-body convergence failed '
|
write(*,*)' Two-body convergence failed '
|
||||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
stop
|
!stop
|
||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
@ -624,13 +624,17 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
|
|||||||
write(*,*) 'Building self-energy'
|
write(*,*) 'Building self-energy'
|
||||||
|
|
||||||
call wall_time(start_t)
|
call wall_time(start_t)
|
||||||
call R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,eOld,EcGM,SigC,Z)
|
call R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eOld,ERI, &
|
||||||
|
eh_sing_rho,old_eh_sing_Om,eh_trip_rho,old_eh_trip_Om, &
|
||||||
|
ee_sing_rho,old_ee_sing_Om,ee_trip_rho,old_ee_trip_Om, &
|
||||||
|
hh_sing_rho,old_hh_sing_Om,hh_trip_rho,old_hh_trip_Om, &
|
||||||
|
EcGM,SigC,Z)
|
||||||
call wall_time(end_t)
|
call wall_time(end_t)
|
||||||
t = end_t - start_t
|
t = end_t - start_t
|
||||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds'
|
write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
eQPlin(:) = eHF(:) !+ Z(:)*SigC(:)
|
eQPlin(:) = eHF(:) + Z(:)*SigC(:)
|
||||||
|
|
||||||
! Solve the quasi-particle equation
|
! Solve the quasi-particle equation
|
||||||
|
|
||||||
@ -647,15 +651,19 @@ subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_i
|
|||||||
write(*,*)
|
write(*,*)
|
||||||
stop
|
stop
|
||||||
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
deallocate(eQPlin,Z,SigC)
|
|
||||||
|
|
||||||
! Check one-body converge
|
! Check one-body converge
|
||||||
|
|
||||||
err_1b = maxval(abs(eOld - eQP))
|
err_1b = maxval(abs(eOld - eQP))
|
||||||
eOld(:) = eQP(:)
|
eOld(:) = eQP(:)
|
||||||
|
|
||||||
|
! Print for one-body part
|
||||||
|
|
||||||
|
call R_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,ERHF,EcGM,Ec_eh,Ec_pp)
|
||||||
|
|
||||||
|
deallocate(eQPlin,Z,SigC)
|
||||||
|
|
||||||
call wall_time(end_1b)
|
call wall_time(end_1b)
|
||||||
t_1b = end_1b - start_1b
|
t_1b = end_1b - start_1b
|
||||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for one-body iteration =',t_1b,' seconds'
|
write(*,'(A50,1X,F9.3,A8)') 'Wall time for one-body iteration =',t_1b,' seconds'
|
||||||
|
|||||||
102
src/Parquet/R_Parquet_self_energy.f90
Normal file
102
src/Parquet/R_Parquet_self_energy.f90
Normal file
@ -0,0 +1,102 @@
|
|||||||
|
subroutine R_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,eQP,ERI, &
|
||||||
|
eh_sing_rho,eh_sing_Om,eh_trip_rho,eh_trip_Om, &
|
||||||
|
ee_sing_rho,ee_sing_Om,ee_trip_rho,ee_trip_Om, &
|
||||||
|
hh_sing_rho,hh_sing_Om,hh_trip_rho,hh_trip_Om, &
|
||||||
|
EcGM,SigC,Z)
|
||||||
|
|
||||||
|
! Compute correlation part of the self-energy with only irreducible vertices contribution
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
integer,intent(in) :: nOrb,nC,nO,nV,nR
|
||||||
|
integer,intent(in) :: nS,nOOs,nVVs,nOOt,nVVt
|
||||||
|
double precision,intent(in) :: eQP(nOrb)
|
||||||
|
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||||
|
double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS)
|
||||||
|
double precision,intent(in) :: eh_sing_Om(nS)
|
||||||
|
double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS)
|
||||||
|
double precision,intent(in) :: eh_trip_Om(nS)
|
||||||
|
double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs)
|
||||||
|
double precision,intent(in) :: ee_sing_Om(nVVs)
|
||||||
|
double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt)
|
||||||
|
double precision,intent(in) :: ee_trip_Om(nVVt)
|
||||||
|
double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs)
|
||||||
|
double precision,intent(in) :: hh_sing_Om(nOOs)
|
||||||
|
double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nOOt)
|
||||||
|
double precision,intent(in) :: hh_trip_Om(nOOt)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
integer :: i,j,k,a,b,c
|
||||||
|
integer :: p,n
|
||||||
|
double precision :: eps,dem1,dem2,reg,reg1,reg2
|
||||||
|
double precision :: num
|
||||||
|
double precision :: start_t,end_t,t
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
double precision,intent(out) :: EcGM
|
||||||
|
double precision,intent(out) :: SigC(nOrb)
|
||||||
|
double precision,intent(out) :: Z(nOrb)
|
||||||
|
|
||||||
|
! Initialize
|
||||||
|
|
||||||
|
SigC(:) = 0d0
|
||||||
|
Z(:) = 0d0
|
||||||
|
EcGM = 0d0
|
||||||
|
|
||||||
|
!-----------------------------!
|
||||||
|
! GF2 part of the self-energy !
|
||||||
|
!-----------------------------!
|
||||||
|
call wall_time(start_t)
|
||||||
|
do p=nC+1,nOrb-nR
|
||||||
|
! 2h1p sum
|
||||||
|
do i=nC+1,nO
|
||||||
|
do j=nC+1,nO
|
||||||
|
do a=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
eps = eQP(p) + eQP(a) - eQP(i) - eQP(j)
|
||||||
|
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
|
||||||
|
num = ERI(p,a,j,i)*(2d0*ERI(j,i,p,a) - ERI(j,i,a,p))
|
||||||
|
|
||||||
|
SigC(p) = SigC(p) + num*reg/eps
|
||||||
|
Z(p) = Z(p) - num*reg/eps**2
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 2p1h sum
|
||||||
|
do i=nC+1,nO
|
||||||
|
do a=nO+1,nOrb-nR
|
||||||
|
do b=nO+1,nOrb-nR
|
||||||
|
|
||||||
|
eps = eQP(p) + eQP(i) - eQP(a) - eQP(b)
|
||||||
|
reg = (1d0 - exp(- 2d0 * eta * eps * eps))
|
||||||
|
num = ERI(p,i,b,a)*(2d0*ERI(b,a,p,i) - ERI(b,a,i,p))
|
||||||
|
|
||||||
|
SigC(p) = SigC(p) + num*reg/eps
|
||||||
|
Z(p) = Z(p) - num*reg/eps**2
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
call wall_time(end_t)
|
||||||
|
t = end_t - start_t
|
||||||
|
|
||||||
|
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
!-----------------------------!
|
||||||
|
! Renormalization factor !
|
||||||
|
!-----------------------------!
|
||||||
|
|
||||||
|
Z(:) = 1d0/(1d0 - Z(:))
|
||||||
|
|
||||||
|
!-------------------------------------!
|
||||||
|
! Galitskii-Migdal correlation energy !
|
||||||
|
!-------------------------------------!
|
||||||
|
|
||||||
|
EcGM = 0d0
|
||||||
|
|
||||||
|
end subroutine
|
||||||
@ -1,94 +0,0 @@
|
|||||||
subroutine R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,e,EcGM,SigC,Z)
|
|
||||||
|
|
||||||
! Compute correlation part of the self-energy with only irreducible vertices contribution
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
integer,intent(in) :: nOrb,nC,nO,nV,nR
|
|
||||||
double precision,intent(in) :: e(nOrb)
|
|
||||||
! Local variables
|
|
||||||
integer :: p,i,j,a,b
|
|
||||||
double precision :: D2p1h,D2h1p
|
|
||||||
! Output variables
|
|
||||||
double precision,intent(out) :: EcGM
|
|
||||||
double precision,intent(out) :: SigC(nOrb)
|
|
||||||
double precision,intent(out) :: Z(nOrb)
|
|
||||||
|
|
||||||
!----------------------------!
|
|
||||||
! Static Parquet self-energy !
|
|
||||||
!----------------------------!
|
|
||||||
SigC(:) = 0d0
|
|
||||||
! 2h1p part of the correlation self-energy
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
do i=nC+1,nO
|
|
||||||
do j=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
D2h1p = e(p) + e(a) - e(i) - e(j)
|
|
||||||
SigC(p) = SigC(p) !+ 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2p1h part of the correlation self-energy
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
do i=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
do b=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
D2p1h = e(p) + e(i) - e(a) - e(b)
|
|
||||||
SigC(p) = SigC(p) !+ 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
!------------------------!
|
|
||||||
! Renormalization factor !
|
|
||||||
!------------------------!
|
|
||||||
Z(:) = 0d0
|
|
||||||
! 2h1p part of the renormlization factor
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
do i=nC+1,nO
|
|
||||||
do j=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
D2h1p = e(p) + e(a) - e(i) - e(j)
|
|
||||||
Z(p) = Z(p)
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2p1h part of the renormlization factor
|
|
||||||
do p=nC+1,nOrb-nR
|
|
||||||
do i=nC+1,nO
|
|
||||||
do a=nO+1,nOrb-nR
|
|
||||||
do b=nO+1,nOrb-nR
|
|
||||||
|
|
||||||
D2p1h = e(p) + e(i) - e(a) - e(b)
|
|
||||||
Z(p) = Z(p)
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
Z(:) = 1d0/(1d0 - Z(:))
|
|
||||||
|
|
||||||
!-------------------------------------!
|
|
||||||
! Galitskii-Migdal correlation energy !
|
|
||||||
!-------------------------------------!
|
|
||||||
|
|
||||||
EcGM = 0d0
|
|
||||||
! do i=nC+1,nO
|
|
||||||
! do a=nO+1,nOrb-nR
|
|
||||||
! do m=1,nS
|
|
||||||
|
|
||||||
! end do
|
|
||||||
! end do
|
|
||||||
! end do
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
64
src/Parquet/R_print_parquet_1b.f90
Normal file
64
src/Parquet/R_print_parquet_1b.f90
Normal file
@ -0,0 +1,64 @@
|
|||||||
|
subroutine R_print_parquet_1b(nOrb,nO,eHF,SigC,eQP,Z,n_it_1b,err_1b,ENuc,EGHF,EcGM,Ec_eh,Ec_pp)
|
||||||
|
|
||||||
|
! Print one-electron energies and other stuff for G0F2
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
integer,intent(in) :: nOrb
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
double precision,intent(in) :: eHF(nOrb)
|
||||||
|
double precision,intent(in) :: SigC(nOrb)
|
||||||
|
double precision,intent(in) :: eQP(nOrb)
|
||||||
|
double precision,intent(in) :: Z(nOrb)
|
||||||
|
integer,intent(in) :: n_it_1b
|
||||||
|
double precision,intent(in) :: err_1b
|
||||||
|
double precision,intent(in) :: ENuc
|
||||||
|
double precision,intent(in) :: EGHF
|
||||||
|
double precision,intent(in) :: EcGM
|
||||||
|
double precision,intent(in) :: Ec_eh(nspin)
|
||||||
|
double precision,intent(in) :: Ec_pp(nspin)
|
||||||
|
|
||||||
|
integer :: p
|
||||||
|
integer :: HOMO
|
||||||
|
integer :: LUMO
|
||||||
|
double precision :: Gap
|
||||||
|
|
||||||
|
! HOMO and LUMO
|
||||||
|
|
||||||
|
HOMO = nO
|
||||||
|
LUMO = HOMO + 1
|
||||||
|
Gap = eQP(LUMO) - eQP(HOMO)
|
||||||
|
|
||||||
|
! Dump results
|
||||||
|
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)' Parquet self-energy '
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||||
|
'|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
|
||||||
|
do p=1,nOrb
|
||||||
|
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||||
|
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eQP(p)*HaToeV,'|'
|
||||||
|
end do
|
||||||
|
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A60,I15)') 'One-body iteration # ',n_it_1b
|
||||||
|
write(*,'(2X,A60,F15.6)') 'One-body convergence ',err_1b
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') 'Parquet HOMO energy = ',eQP(HOMO)*HaToeV,' eV'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') 'Parquet LUMO energy = ',eQP(LUMO)*HaToeV,' eV'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') 'Parquet HOMO-LUMO gap = ',Gap*HaToeV,' eV'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') ' Parquet total energy = ',ENuc + EGHF + EcGM,' au'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') ' Parquet correlation energy = ',EcGM,' au'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') ' eh-RPA correlation energy = ',Ec_eh(1)+3d0*Ec_eh(2),' au'
|
||||||
|
write(*,'(2X,A60,F15.6,A3)') ' pp-RPA correlation energy = ',Ec_pp(1)+3d0*Ec_pp(2),' au'
|
||||||
|
!write(*,'(2X,A60,F15.6,A3)') '(eh+pp)-RPA correlation energy = ',Ec_pp,' au'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end subroutine
|
||||||
Loading…
x
Reference in New Issue
Block a user