LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

cleaning up ufGW

This commit is contained in:
Pierre-Francois Loos 2023-11-24 15:31:29 +01:00
parent 3f12eb8188
commit 5d74b8073a
4 changed files with 595 additions and 404 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/GW/GGW.f90
View File

@ -1,6 +1,6 @@
subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
TDA_W,TDA,dBSE,dTDA,linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
TDA_W,TDA,dBSE,dTDA,linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! GW module
@ -45,8 +45,8 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas2)
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas2)
double precision,intent(in) :: cHF(nBas2,nBas2)
double precision,intent(in) :: PHF(nBas2,nBas2)
double precision,intent(in) :: S(nBas2,nBas2)
@ -71,7 +71,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
call wall_time(start_GW)
call GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -88,7 +88,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
call wall_time(start_GW)
call evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -105,8 +105,8 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
call wall_time(start_GW)
call qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW

28
src/GW/RGW.f90
View File

@ -1,7 +1,7 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! GW module
@ -50,8 +50,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
@ -76,7 +76,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -93,7 +93,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -110,8 +110,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -128,8 +128,8 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -145,7 +145,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufG0W0) then
call wall_time(start_GW)
call ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF,TDA_W)
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -161,7 +161,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufGW) then
call wall_time(start_GW)
call ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF)
call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW

538
src/GW/ufG0W0.f90
View File

@ -1,4 +1,4 @@
subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold G0W0 equations
@ -9,6 +9,7 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
logical,intent(in) :: dotest
logical,intent(in) :: TDA_W
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -19,7 +20,6 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
logical,intent(in) :: TDA_W
! Local variables
@ -36,7 +36,6 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
integer :: n2h1p,n2p1h,nH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
double precision,allocatable :: cGW(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: Aph(:,:)
@ -55,9 +54,9 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
! Hello world
write(*,*)
write(*,*)'**********************************************'
write(*,*)'| Unfolded G0W0 calculation |'
write(*,*)'**********************************************'
write(*,*)'****************************************'
write(*,*)'* Restricted Upfolded G0W0 Calculation *'
write(*,*)'****************************************'
write(*,*)
! Dimension of the supermatrix
@ -68,7 +67,7 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
! Memory allocation
allocate(H(nH,nH),cGW(nH,nH),eGW(nH),Z(nH))
allocate(H(nH,nH),eGW(nH),Z(nH))
! Initialization
@ -77,230 +76,230 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
H(:,:) = 0d0
!!! Compute only the HOMO !!!
p=nO
p = nO
if (TDA_W) then
! TDA for W
! TDA for W
write(*,*) 'Tamm-Dancoff approximation actived!'
write(*,*)
write(*,*) 'Tamm-Dancoff approximation actived!'
write(*,*)
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
!-------------!
! Block C2h1p !
!-------------!
!---------!
! Block F !
!---------!
H(1,1) = eHF(p)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
H(1+ija,1+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)
end do
end do
end do
!-------------!
! Block V2h1p !
!-------------!
end do
end do
end do
!-------------!
! Block C2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
H(1+n2h1p+iab,1+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)
end do
end do
end do
end do
end do
end do
!---------!
! Block F !
!---------!
H(1,1) = eHF(p)
!-------------!
! Block V2h1p !
!-------------!
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
H(1 ,1+klc) = sqrt(2d0)*ERI(p,c,k,l)
H(1+klc,1 ) = sqrt(2d0)*ERI(p,c,k,l)
end do
end do
end do
!-------------!
! Block V2p1h !
!-------------!
kcd = 0
do k=nC+1,nO
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
H(1 ,1+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
H(1+n2h1p+kcd,1 ) = sqrt(2d0)*ERI(p,k,d,c)
end do
klc = klc + 1
H(1 ,1+klc) = sqrt(2d0)*ERI(p,c,k,l)
H(1+klc,1 ) = sqrt(2d0)*ERI(p,c,k,l)
end do
end do
end do
end do
!-------------!
! Block V2p1h !
!-------------!
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
H(1 ,1+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
H(1+n2h1p+kcd,1 ) = sqrt(2d0)*ERI(p,k,d,c)
end do
end do
end do
!-------------!
! Block C2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
H(1+ija,1+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)
end do
end do
end do
end do
end do
end do
!-------------!
! Block C2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
H(1+n2h1p+iab,1+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)
end do
end do
end do
end do
end do
end do
else
! RPA for W
! RPA for W
write(*,*) 'Tamm-Dancoff approximation deactivated!'
write(*,*)
write(*,*) 'Tamm-Dancoff approximation deactivated!'
write(*,*)
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F W2h1p W2p1h | !
! | | !
! H = | W2h1p D2h1p 0 | !
! | | !
! | W2p1h 0 D2p1h | !
! !
!---------------------------!
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F W2h1p W2p1h | !
! | | !
! H = | W2h1p D2h1p 0 | !
! | | !
! | W2p1h 0 D2p1h | !
! !
!---------------------------!
! Memory allocation !
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS))
! Memory allocation
! Spin manifold
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS))
ispin = 1
! Spin manifold
!-------------------!
! Compute screening !
!-------------------!
ispin = 1
!-------------------!
! Compute screening !
!-------------------!
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------------------!
! Compute spectral weights !
!--------------------------!
!--------------------------!
! Compute spectral weights !
!--------------------------!
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
!---------!
! Block F !
!---------!
!---------!
! Block F !
!---------!
H(1,1) = eHF(p)
H(1,1) = eHF(p)
!-------------!
! Block D2h1p !
!-------------!
!-------------!
! Block D2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
H(1+ija,1+ija) = eHF(i) - Om(ja)
end do
end do
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
!-------------!
! Block W2h1p !
!-------------!
H(1+ija,1+ija) = eHF(i) - Om(ja)
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja)
H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja)
end do
end do
end do
end do
!-------------!
! Block D2h1p !
!-------------!
!-------------!
! Block W2h1p !
!-------------!
iab = 0
do b=nO+1,nBas-nR
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
iab = iab + 1
H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia)
end do
end do
end do
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
!-------------!
! Block W2p1h !
!-------------!
H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja)
H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja)
iab = 0
do b=nO+1,nBas-nR
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
iab = iab + 1
H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia)
H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia)
end do
end do
end do
end do
end do
!-------------!
! Block D2p1h !
!-------------!
iab = 0
do b=nO+1,nBas-nR
do ia=1,nS
iab = iab + 1
H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia)
end do
end do
!-------------!
! Block W2p1h !
!-------------!
iab = 0
do b=nO+1,nBas-nR
do ia=1,nS
iab = iab + 1
H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia)
H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia)
end do
end do
end if
@ -308,96 +307,95 @@ subroutine ufG0W0(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
! Diagonalize supermatrix !
!-------------------------!
cGW(:,:) = H(:,:)
call diagonalize_matrix(nH,cGW,eGW)
H(:,:) = H(:,:)
call diagonalize_matrix(nH,H,eGW)
!-----------------!
! Compute weights !
!-----------------!
do s=1,nH
Z(s) = cGW(1,s)**2
Z(s) = H(1,s)**2
end do
!--------------!
! Dump results !
!--------------!
write(*,*)'-------------------------------------------'
write(*,'(A35,I3)')' G0W0 energies (eV) for orbital ',p
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP (eV)','|','Z','|'
write(*,*)'-------------------------------------------'
write(*,*)'-------------------------------------------'
write(*,'(1X,A32,I3,A8)')'| G0W0 energies (eV) for orbital',p,' |'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP','|','Z','|'
write(*,*)'-------------------------------------------'
do s=1,nH
if(Z(s) > cutoff1) then
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
end if
end do
write(*,*)'-------------------------------------------'
write(*,*)
if(verbose) then
do s=1,nH
if(Z(s) > cutoff1) then
write(*,*)'*************************************************************'
write(*,'(1X,A20,I3,A6,I3)')'Vector for orbital ',p,' and #',s
write(*,'(1X,A7,F10.6,A13,F10.6,1X)')' e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s)
write(*,*)'*************************************************************'
write(*,'(1X,A20,1X,A20,1X,A15,1X)') &
' Configuration ',' Coefficient ',' Weight '
write(*,*)'*************************************************************'
if(p <= nO) &
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,') ',cGW(1,s),cGW(1,s)**2
if(p > nO) &
write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') &
' (',p,') ',cGW(1,s),cGW(1,s)**2
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
if(abs(cGW(1+klc,s)) > cutoff2) &
write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') &
' (',k,',',l,') -> (',c,') ',cGW(1+klc,s),cGW(1+klc,s)**2
end do
end do
end do
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
if(abs(cGW(1+n2h1p+kcd,s)) > cutoff2) &
write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') &
' (',k,') -> (',c,',',d,') ',cGW(1+n2h1p+kcd,s),cGW(1+n2h1p+kcd,s)**2
end do
end do
end do
write(*,*)'*************************************************************'
write(*,*)
end if
end do
do s=1,nH
if(Z(s) > cutoff1) then
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
end if
end do
write(*,*)'-------------------------------------------'
write(*,*)
if(verbose) then
do s=1,nH
if(Z(s) > cutoff1) then
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') &
'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A20,1X,A20,1X,A15,1X)') &
' Configuration ',' Coefficient ',' Weight '
write(*,*)'-------------------------------------------------------------'
if(p <= nO) &
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(1,s),H(1,s)**2
if(p > nO) &
write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(1,s),H(1,s)**2
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
if(abs(H(1+klc,s)) > cutoff2) &
write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') &
' (',k,',',l,') -> (',c,') ',H(1+klc,s),H(1+klc,s)**2
end do
end do
end do
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
if(abs(H(1+n2h1p+kcd,s)) > cutoff2) &
write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') &
' (',k,') -> (',c,',',d,') ',H(1+n2h1p+kcd,s),H(1+n2h1p+kcd,s)**2
end do
end do
end do
write(*,*)'-------------------------------------------------------------'
write(*,*)
end if
end do
end if
end subroutine

415
src/GW/ufGW.f90
View File

@ -1,4 +1,4 @@
subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold GW equations
@ -9,6 +9,7 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
logical,intent(in) :: dotest
logical,intent(in) :: TDA_W
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -26,27 +27,36 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
integer :: s
integer :: i,j,k,l
integer :: a,b,c,d
integer :: ia,ja,kc,lc
integer :: klc,kcd,ija,iab
logical :: dRPA
integer :: ispin
double precision :: EcRPA
integer :: n2h1p,n2p1h,nH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
logical :: verbose = .true.
double precision,parameter :: cutoff1 = 0.01d0
double precision,parameter :: cutoff2 = 0.01d0
! 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(*,*)'**************************************'
write(*,*)'* Restricted Upfolded GW Calculation *'
write(*,*)'**************************************'
write(*,*)
! Dimension of the supermatrix
@ -63,123 +73,248 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
H(:,:) = 0d0
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
if (TDA_W) then
!---------!
! Block F !
!---------!
! TDA for W
do p=nC+1,nBas-nR
H(p,p) = eHF(p)
end do
write(*,*) 'Tamm-Dancoff approximation actived!'
write(*,*)
!-------------!
! Block V2h1p !
!-------------!
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
do p=nC+1,nBas-nR
klc = 0
do k=nC+1,nO
do l=nC+1,nO
!---------!
! Block F !
!---------!
do p=nC+1,nBas-nR
H(p,p) = eHF(p)
end do
!-------------!
! Block V2h1p !
!-------------!
do p=nC+1,nBas-nR
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
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 !
!-------------!
do p=nC+1,nBas-nR
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*ERI(p,k,d,c)
end do
end do
end do
end do
!-------------!
! Block C2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
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)
end do
end do
end do
end do
end do
end do
!-------------!
! Block C2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
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)
end do
end do
end do
end do
end do
end do
else
! RPA for W
write(*,*) 'Tamm-Dancoff approximation deactivated!'
write(*,*)
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F W2h1p W2p1h | !
! | | !
! H = | W2h1p D2h1p 0 | !
! | | !
! | W2p1h 0 D2p1h | !
! !
!---------------------------!
! Memory allocation
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS))
! Spin manifold
ispin = 1
!-------------------!
! Compute screening !
!-------------------!
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------------------!
! Compute spectral weights !
!--------------------------!
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
!---------!
! Block F !
!---------!
do p=nC+1,nBas-nR
H(p,p) = eHF(p)
end do
!-------------!
! Block W2h1p !
!-------------!
do p=nC+1,nBas-nR
klc = 0
do k=nC+1,nO
do lc=1,nS
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)
H(p ,nBas+klc) = sqrt(2d0)*rho(p,k,lc)
H(nBas+klc,p ) = sqrt(2d0)*rho(p,k,lc)
end do
end do
end do
end do
!-------------!
! Block V2p1h !
!-------------!
do p=nC+1,nBas-nR
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
!-------------!
! Block W2p1h !
!-------------!
do p=nC+1,nBas-nR
kcd = 0
do kc=1,nS
do d=nO+1,nBas-nR
kcd = kcd + 1
H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*ERI(p,k,d,c)
H(p ,nBas+n2h1p+kcd) = sqrt(2d0)*rho(p,d,kc)
H(nBas+n2h1p+kcd,p ) = sqrt(2d0)*rho(p,d,kc)
end do
end do
end do
end do
!-------------!
! Block C2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
!-------------!
! Block D2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
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)
end do
end do
end do
H(nBas+ija,nBas+ija) = eHF(i) - Om(ja)
end do
end do
end do
!-------------!
! Block C2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
!-------------!
! Block D2p1h !
!-------------!
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
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)
end do
end do
end do
H(nBas+n2h1p+iab,nBas+n2h1p+iab) = eHF(b) + Om(ia)
end do
end do
end do
end if
!-------------------------!
! Diagonalize supermatrix !
@ -202,19 +337,77 @@ subroutine ufGW(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Dump results !
!--------------!
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(*,*)'-------------------------------------------'
write(*,*)'---------------------------------------------'
write(*,'(1X,A45)')'| GW energies (eV) for all orbitals |'
write(*,*)'---------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP','|','Z','|'
write(*,*)'---------------------------------------------'
do s=1,nH
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
write(*,'(1X,A1,1X,I5,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
enddo
write(*,*)'-------------------------------------------'
write(*,*)'---------------------------------------------'
write(*,*)
if(verbose) then
do s=1,nH
if(Z(s) > cutoff1) then
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A10,I5,A1,1X,A7,F12.6,A13,F6.4,1X)') &
'Solution',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A20,1X,A20,1X,A15,1X)') &
' Configuration ',' Coefficient ',' Weight '
write(*,*)'-------------------------------------------------------------'
if(p <= nO) &
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(1,s),H(1,s)**2
if(p > nO) &
write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(1,s),H(1,s)**2
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
klc = klc + 1
if(abs(H(1+klc,s)) > cutoff2) &
write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') &
' (',k,',',l,') -> (',c,') ',H(1+klc,s),H(1+klc,s)**2
end do
end do
end do
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
kcd = kcd + 1
if(abs(H(1+n2h1p+kcd,s)) > cutoff2) &
write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') &
' (',k,') -> (',c,',',d,') ',H(1+n2h1p+kcd,s),H(1+n2h1p+kcd,s)**2
end do
end do
end do
write(*,*)'-------------------------------------------------------------'
write(*,*)
end if
end do
end if
end subroutine