10
1
mirror of https://github.com/pfloos/quack synced 2024-12-25 13:53:41 +01:00
QuAcK/src/GW/ufGW.f90

221 lines
4.6 KiB
Fortran
Raw Normal View History

2021-12-17 11:41:40 +01:00
subroutine ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
2021-10-18 10:10:54 +02:00
! Unfold GW equations
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
! Local variables
integer :: p
2021-10-18 15:05:01 +02:00
integer :: s
2021-10-18 10:10:54 +02:00
integer :: i,j,k,l
integer :: a,b,c,d
integer :: klc,kcd,ija,iab
integer :: n2h1p,n2p1h,nH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
double precision,allocatable :: eGW(:)
2021-10-18 15:05:01 +02:00
double precision,allocatable :: Z(:)
2021-10-18 10:10:54 +02:00
! Output variables
! Hello world
write(*,*)
write(*,*)'**********************************************'
write(*,*)'| Unfolded GW calculation |'
write(*,*)'**********************************************'
write(*,*)
! TDA for W
write(*,*) 'Tamm-Dancoff approximation for dynamic screening by default!'
write(*,*)
! Dimension of the supermatrix
2023-04-12 21:41:30 +02:00
n2h1p = nO*nO*nV
2021-10-18 14:14:46 +02:00
n2p1h = nV*nV*nO
2021-10-18 10:10:54 +02:00
nH = nBas + n2h1p + n2p1h
! Memory allocation
2021-10-18 15:05:01 +02:00
allocate(H(nH,nH),eGW(nH),Z(nH))
2021-10-18 10:10:54 +02:00
! Initialization
H(:,:) = 0d0
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
!---------!
! Block F !
!---------!
2021-10-18 14:14:46 +02:00
do p=nC+1,nBas-nR
2021-10-18 10:10:54 +02:00
H(p,p) = eHF(p)
end do
!-------------!
! Block V2h1p !
!-------------!
2021-10-18 14:14:46 +02:00
do p=nC+1,nBas-nR
2021-10-18 10:10:54 +02:00
klc = 0
do k=nC+1,nO
do l=nC+1,nO
2021-10-18 14:14:46 +02:00
do c=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
klc = klc + 1
H(p ,nBas+klc) = sqrt(2d0)*ERI(p,c,k,l)
H(nBas+klc,p ) = sqrt(2d0)*ERI(p,c,k,l)
end do
end do
end do
end do
!-------------!
! Block V2p1h !
!-------------!
2021-10-18 14:14:46 +02:00
do p=nC+1,nBas-nR
2021-10-18 10:10:54 +02:00
kcd = 0
do k=nC+1,nO
2021-10-18 14:14:46 +02:00
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
kcd = kcd + 1
2021-10-18 21:11:01 +02:00
H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
2021-10-18 14:14:46 +02:00
H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*ERI(p,k,d,c)
2021-10-18 10:10:54 +02:00
end do
end do
end do
end do
!-------------!
! Block C2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
2021-10-18 14:14:46 +02:00
do a=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
2021-10-18 14:14:46 +02:00
do c=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
klc = klc + 1
2021-10-18 14:14:46 +02:00
H(nBas+ija,nBas+klc) &
= ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
- 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k)
2021-10-18 10:10:54 +02:00
end do
end do
end do
end do
end do
end do
!-------------!
! Block C2p1h !
!-------------!
iab = 0
do i=nC+1,nO
2021-10-18 14:14:46 +02:00
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
iab = iab + 1
kcd = 0
do k=nC+1,nO
2021-10-18 14:14:46 +02:00
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
2021-10-18 10:10:54 +02:00
kcd = kcd + 1
2021-10-18 14:14:46 +02:00
H(nBas+n2h1p+iab,nBas+n2h1p+kcd) &
= ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
+ 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d)
2021-10-18 10:10:54 +02:00
end do
end do
end do
end do
end do
end do
!-------------------------!
! Diagonalize supermatrix !
!-------------------------!
call diagonalize_matrix(nH,H,eGW)
2021-10-18 15:05:01 +02:00
!-----------------!
! Compute weights !
!-----------------!
Z(:) = 0d0
do s=1,nH
do p=nC+1,nBas-nR
Z(s) = Z(s) + H(p,s)**2
end do
end do
2021-10-18 10:10:54 +02:00
!--------------!
! Dump results !
!--------------!
2021-10-18 21:11:01 +02:00
write(*,*)'-------------------------------------------'
write(*,*)' unfolded GW energies (eV) '
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP (eV)','|','Z','|'
write(*,*)'-------------------------------------------'
do s=1,nH
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
enddo
write(*,*)'-------------------------------------------'
write(*,*)
2021-10-18 10:10:54 +02:00
end subroutine ufGW