LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-09 07:33:55 +01:00
Code Releases Activity

moving calculations of Z in Sig

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-12 14:13:45 +02:00
parent aec99cb41e
commit d6eb90df47
18 changed files with 263 additions and 371 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/GT/G0T0eh.f90
View File

@ -53,7 +53,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
double precision :: EcppBSE(nspin) double precision :: EcppBSE(nspin)
double precision :: EcGM double precision :: EcGM
double precision,allocatable :: SigX(:) double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
@ -99,7 +99,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
! Memory allocation ! Memory allocation
allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & allocate(Sig(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas)) rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas))
!-------------------! !-------------------!
@ -124,13 +124,12 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
if(regularize) then if(regularize) then
! call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) ! call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,Sig)
! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) ! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
else else
call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z)
call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,Z)
end if end if
@ -138,7 +137,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
! Solve the quasi-particle equation ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:))
! Linearized or graphical solution? ! Linearized or graphical solution?
@ -170,11 +169,11 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
! Dump results ! ! Dump results !
!--------------! !--------------!
call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
! Deallocate memory ! Deallocate memory
! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) ! deallocate(Sig,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin)
! Plot stuff ! Plot stuff

20
src/GT/G0T0pp.f90
View File

@ -60,7 +60,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
double precision,allocatable :: Y2ab(:,:),Y2aa(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:)
double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:)
double precision,allocatable :: SigX(:) double precision,allocatable :: SigX(:)
double precision,allocatable :: SigT(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: eGTlin(:) double precision,allocatable :: eGTlin(:)
@ -92,7 +92,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), & Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), &
Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), &
rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), &
SigX(nBas),SigT(nBas),Z(nBas),eGTlin(nBas)) SigX(nBas),Sig(nBas),Z(nBas),eGTlin(nBas))
!---------------------------------------------- !----------------------------------------------
! alpha-beta block ! alpha-beta block
@ -127,7 +127,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
!---------------------------------------------- !----------------------------------------------
EcGM = 0d0 EcGM = 0d0
SigT(:) = 0d0 Sig(:) = 0d0
Z(:) = 0d0 Z(:) = 0d0
iblock = 3 iblock = 3
@ -136,13 +136,12 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
if(regularize) then if(regularize) then
call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,SigT) call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,Sig)
call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,Z) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,Z)
else else
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,SigT) call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,Sig,Z)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,Z)
end if end if
@ -152,13 +151,12 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
if(regularize) then if(regularize) then
call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,SigT) call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,Sig)
call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,Z) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,Z)
else else
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,SigT) call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,Sig,Z)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,Z)
end if end if
@ -174,7 +172,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
! Solve the quasi-particle equation ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + SigT(:) - Vxc(:)) eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:))
if(linearize) then if(linearize) then
@ -210,7 +208,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2) EcRPA(2) = 3d0*EcRPA(2)
call print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) call print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
! Perform BSE calculation ! Perform BSE calculation

61
src/GT/GTeh_renormalization_factor.f90
View File

@ -1,61 +0,0 @@
subroutine GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,Z)
! Compute renormalization factor for GW
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
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) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS)
double precision,intent(in) :: rhoR(nBas,nBas,nS)
! Local variables
integer :: p,i,a,m
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do m=1,nS
eps = e(p) - e(i) + Om(m)
Z(p) = Z(p) - rhoL(i,p,m)*rhoR(i,p,m)*(eps/(eps**2 + eta**2))**2
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do m=1,nS
eps = e(p) - e(a) - Om(m)
Z(p) = Z(p) - rhoL(p,a,m)*rhoR(p,a,m)*(eps/(eps**2 + eta**2))**2
end do
end do
end do
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine

39
src/GT/GTeh_self_energy.f90
View File

@ -1,6 +1,6 @@
subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig,Z)
! Compute correlation part of the self-energy for GTeh ! Compute correlation part of the self-energy for GTeh and the renormalization factor
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -24,16 +24,18 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC)
integer :: i,j,a,b integer :: i,j,a,b
integer :: p,q,r integer :: p,q,r
integer :: m integer :: m
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(out) :: EcGM double precision,intent(out) :: EcGM
double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: Sig(nBas,nBas)
double precision,intent(out) :: Z(nBas)
! Initialize ! Initialize
SigC(:,:) = 0d0 Sig(:,:) = 0d0
Z(:) = 0d0
!----------------! !----------------!
! GW self-energy ! ! GW self-energy !
@ -42,16 +44,20 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC)
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(m,i,q,p,eps) & !$OMP PRIVATE(m,i,q,p,num,eps) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do m=1,nS do m=1,nS
do i=nC+1,nO do i=nC+1,nO
eps = e(p) - e(i) + Om(m) eps = e(p) - e(i) + Om(m)
SigC(p,q) = SigC(p,q) + rhoL(i,p,m)*rhoR(i,q,m)*eps/(eps**2 + eta**2) num = rhoL(i,p,m)*rhoR(i,q,m)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -62,16 +68,20 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC)
! Virtual part of the correlation self-energy ! Virtual part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(m,a,q,p,eps) & !$OMP PRIVATE(m,a,q,p,num,eps) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do m=1,nS do m=1,nS
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
eps = e(p) - e(a) - Om(m) eps = e(p) - e(a) - Om(m)
SigC(p,q) = SigC(p,q) + rhoL(p,a,m)*rhoR(q,a,m)*eps/(eps**2 + eta**2) num = rhoL(p,a,m)*rhoR(q,a,m)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -86,9 +96,14 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC)
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
eps = e(a) - e(i) + Om(m) eps = e(a) - e(i) + Om(m)
EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2) num = rhoL(i,a,m)*rhoR(i,a,m)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
end do end do
end do end do
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

33
src/GT/GTeh_self_energy_diag.f90
View File

@ -1,6 +1,6 @@
subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the self-energy ! Compute diagonal of the correlation part of the self-energy and the renormalization factor
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -22,16 +22,18 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
! Local variables ! Local variables
integer :: i,a,p,q,m integer :: i,a,p,q,m
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(out) :: SigC(nBas) double precision,intent(out) :: Sig(nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: EcGM double precision,intent(out) :: EcGM
! Initialize ! Initialize
SigC(:) = 0d0 Sig(:) = 0d0
Z(:) = 0d0
!----------------------------- !-----------------------------
! GW self-energy ! GW self-energy
@ -42,8 +44,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do m=1,nS do m=1,nS
eps = e(p) - e(i) + Om(m) eps = e(p) - e(i) + Om(m)
SigC(p) = SigC(p) + rhoL(i,p,m)*rhoR(i,p,m)*eps/(eps**2 + eta**2) num = rhoL(i,p,m)*rhoR(i,p,m)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -53,8 +59,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do m=1,nS do m=1,nS
eps = e(p) - e(a) - Om(m) eps = e(p) - e(a) - Om(m)
SigC(p) = SigC(p) + rhoL(p,a,m)*rhoR(p,a,m)*eps/(eps**2 + eta**2) num = rhoL(p,a,m)*rhoR(p,a,m)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -65,10 +75,17 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
do i=nC+1,nO do i=nC+1,nO
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do m=1,nS do m=1,nS
eps = e(a) - e(i) + Om(m) eps = e(a) - e(i) + Om(m)
EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2) num = rhoL(i,a,m)*rhoR(i,a,m)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
end do end do
end do end do
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

51
src/GT/GTpp_renormalization_factor.f90
View File

@ -1,51 +0,0 @@
subroutine GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z)
! Compute renormalization factor of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,a,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) - rho1(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=1,nV-nR
do kl=1,nOO
eps = e(p) + e(nO+a) - Omega2(kl)
Z(p) = Z(p) - rho2(p,nO+a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
end subroutine

43
src/GT/GTpp_self_energy.f90
View File

@ -1,6 +1,6 @@
subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
! Compute the correlation part of the T-matrix self-energy ! Compute the correlation part of the T-matrix self-energy and the renormalization factor
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -16,20 +16,21 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO) double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables ! Local variables
integer :: i,j,a,b,p,q,cd,kl integer :: i,j,a,b,p,q,cd,kl
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(inout):: EcGM double precision,intent(inout):: EcGM
double precision,intent(inout):: SigT(nBas,nBas) double precision,intent(inout):: Sig(nBas,nBas)
double precision,intent(inout):: Z(nBas)
!---------------------------------------------- !----------------------------------------------
! Occupied part of the T-matrix self-energy ! Occupied part of the T-matrix self-energy
@ -39,8 +40,12 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do cd=1,nVV do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2) eps = e(p) + e(i) - Om1(cd)
num = rho1(p,i,cd)*rho1(q,i,cd)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo
enddo enddo
@ -54,8 +59,12 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do kl=1,nOO do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2) eps = e(p) + e(a) - Om2(kl)
num = rho2(p,a,kl)*rho2(q,a,kl)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo
enddo enddo
@ -68,8 +77,11 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh
do i=nC+1,nO do i=nC+1,nO
do j=nC+1,nO do j=nC+1,nO
do cd=1,nVV do cd=1,nVV
eps = e(i) + e(j) - Omega1(cd)
EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) eps = e(i) + e(j) - Om1(cd)
num = rho1(i,j,cd)*rho1(i,j,cd)
EcGM = EcGM + num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo
@ -77,8 +89,11 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
do kl=1,nOO do kl=1,nOO
eps = e(a) + e(b) - Omega2(kl)
EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) eps = e(a) + e(b) - Om2(kl)
num = rho2(a,b,kl)*rho2(a,b,kl)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo

41
src/GT/GTpp_self_energy_diag.f90
View File

@ -1,4 +1,4 @@
subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the T-matrix self-energy ! Compute diagonal of the correlation part of the T-matrix self-energy
@ -16,20 +16,21 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO) double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables ! Local variables
integer :: i,j,a,b,p,cd,kl integer :: i,j,a,b,p,cd,kl
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(inout) :: EcGM double precision,intent(inout) :: EcGM
double precision,intent(inout) :: SigT(nBas) double precision,intent(inout) :: Sig(nBas)
double precision,intent(inout) :: Z(nBas)
!---------------------------------------------- !----------------------------------------------
! Occupied part of the T-matrix self-energy ! Occupied part of the T-matrix self-energy
@ -38,8 +39,12 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do cd=1,nVV do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) eps = e(p) + e(i) - Om1(cd)
num = rho1(p,i,cd)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo
enddo enddo
@ -51,8 +56,12 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do kl=1,nOO do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) eps = e(p) + e(a) - Om2(kl)
num = rho2(p,a,kl)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo
enddo enddo
enddo enddo
@ -64,8 +73,11 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg
do i=nC+1,nO do i=nC+1,nO
do j=nC+1,nO do j=nC+1,nO
do cd=1,nVV do cd=1,nVV
eps = e(i) + e(j) - Omega1(cd)
EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) eps = e(i) + e(j) - Om1(cd)
num = rho1(i,j,cd)**2
EcGM = EcGM + num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo
@ -73,8 +85,11 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
do kl=1,nOO do kl=1,nOO
eps = e(a) + e(b) - Omega2(kl)
EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) eps = e(a) + e(b) - Om2(kl)
num = rho2(a,b,kl)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo

15
src/GT/evGTeh.f90
View File

@ -68,7 +68,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
double precision,allocatable :: eOld(:) double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: SigX(:) double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:)
@ -106,7 +106,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Memory allocation ! Memory allocation
allocate(eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & allocate(eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),Sig(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas)) rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas))
! Compute the exchange part of the self-energy ! Compute the exchange part of the self-energy
@ -146,19 +146,18 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
if(regularize) then if(regularize) then
! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) ! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,Sig)
! call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) ! call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
else else
call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z)
call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z)
end if end if
! Solve the quasi-particle equation ! Solve the quasi-particle equation
eGTlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) eGTlin(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:)
! Linearized or graphical solution? ! Linearized or graphical solution?
@ -184,7 +183,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Print results ! Print results
call print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) call print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
! Linear mixing or DIIS extrapolation ! Linear mixing or DIIS extrapolation
@ -232,7 +231,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Deallocate memory ! Deallocate memory
deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error_diis,e_diis) deallocate(eOld,Z,Sig,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error_diis,e_diis)
! Perform BSE calculation ! Perform BSE calculation

28
src/GT/evGTpp.f90
View File

@ -71,7 +71,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
double precision,allocatable :: SigX(:) double precision,allocatable :: SigX(:)
double precision,allocatable :: SigT(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
! Output variables ! Output variables
@ -100,7 +100,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), &
eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigT(nBas), & eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),Sig(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis)) error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Compute the exchange part of the self-energy ! Compute the exchange part of the self-energy
@ -155,30 +155,20 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
!---------------------------------------------- !----------------------------------------------
EcGM = 0d0 EcGM = 0d0
SigT(:) = 0d0 Sig(:) = 0d0
Z(:) = 0d0 Z(:) = 0d0
iblock = 3 iblock = 3
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, & call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig,Z)
Om1s,rho1s,Om2s,rho2s,EcGM,SigT)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
Om1s,rho1s,Om2s,rho2s,Z)
iblock = 4 iblock = 4
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, & call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
Om1t,rho1t,Om2t,rho2t,EcGM,SigT)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
Om1t,rho1t,Om2t,rho2t,Z)
Z(:) = 1d0/(1d0 - Z(:)) Z(:) = 1d0/(1d0 - Z(:))
@ -188,7 +178,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
! Solve the quasi-particle equation ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
eGT(:) = eHF(:) + SigX(:) + SigT(:) - Vxc(:) eGT(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:)
! Convergence criteria ! Convergence criteria
@ -198,7 +188,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
! Dump results ! Dump results
!---------------------------------------------- !----------------------------------------------
call print_evGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) call print_evGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
! DIIS extrapolation ! DIIS extrapolation

27
src/GT/qsGTeh.f90
View File

@ -89,9 +89,9 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
double precision,allocatable :: J(:,:) double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:) double precision,allocatable :: K(:,:)
double precision,allocatable :: SigC(:,:) double precision,allocatable :: Sig(:,:)
double precision,allocatable :: SigCp(:,:) double precision,allocatable :: Sigp(:,:)
double precision,allocatable :: SigCm(:,:) double precision,allocatable :: Sigm(:,:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:) double precision,allocatable :: error(:,:)
@ -129,7 +129,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Memory allocation ! Memory allocation
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), &
OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS), & OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
@ -180,26 +180,25 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
if(regularize) then if(regularize) then
! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) ! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig)
! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) ! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z)
else else
call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z)
call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z)
endif endif
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
SigCp = 0.5d0*(SigC + transpose(SigC)) Sigp = 0.5d0*(Sig + transpose(Sig))
SigCm = 0.5d0*(SigC - transpose(SigC)) Sigm = 0.5d0*(Sig - transpose(Sig))
call MOtoAO_transform(nBas,S,c,SigCp) call MOtoAO_transform(nBas,S,c,Sigp)
! Solve the quasi-particle equation ! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:)
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
@ -220,7 +219,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT) call diagonalize_matrix(nBas,cp,eGT)
c = matmul(X,cp) c = matmul(X,cp)
SigCp = matmul(transpose(c),matmul(SigCp,c)) Sigp = matmul(transpose(c),matmul(Sigp,c))
! Compute new density matrix in the AO basis ! Compute new density matrix in the AO basis
@ -258,7 +257,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Print results ! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
enddo enddo
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -281,7 +280,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis) deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis)
! Perform BSE calculation ! Perform BSE calculation

36
src/GT/qsGTpp.f90
View File

@ -89,9 +89,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
double precision,allocatable :: J(:,:) double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:) double precision,allocatable :: K(:,:)
double precision,allocatable :: SigT(:,:) double precision,allocatable :: Sig(:,:)
double precision,allocatable :: SigTp(:,:) double precision,allocatable :: Sigp(:,:)
double precision,allocatable :: SigTm(:,:) double precision,allocatable :: Sigm(:,:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:) double precision,allocatable :: error(:,:)
@ -137,7 +137,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Memory allocation ! Memory allocation
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigT(nBas,nBas),SigTp(nBas,nBas),SigTm(nBas,nBas),Z(nBas), & J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
@ -201,7 +201,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
EcGM = 0d0 EcGM = 0d0
SigT(:,:) = 0d0 Sig(:,:) = 0d0
Z(:) = 0d0 Z(:) = 0d0
iblock = 3 iblock = 3
@ -210,15 +210,13 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
if(regularize) then if(regularize) then
call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,SigT) call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig)
call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,Z) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,Z)
else else
call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,SigT) call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig,Z)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,Z)
end if end if
@ -228,16 +226,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
if(regularize) then if(regularize) then
call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,SigT) call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig)
call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
Om1t,rho1t,Om2t,rho2t,Z) Om1t,rho1t,Om2t,rho2t,Z)
else else
call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,SigT) call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,Z)
end if end if
@ -245,14 +241,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
SigTp = 0.5d0*(SigT + transpose(SigT)) Sigp = 0.5d0*(Sig + transpose(Sig))
SigTm = 0.5d0*(SigT - transpose(SigT)) Sigm = 0.5d0*(Sig - transpose(Sig))
call MOtoAO_transform(nBas,S,c,SigTp) call MOtoAO_transform(nBas,S,c,Sigp)
! Solve the quasi-particle equation ! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigTp(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:)
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
@ -273,7 +269,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT) call diagonalize_matrix(nBas,cp,eGT)
c = matmul(X,cp) c = matmul(X,cp)
SigTp = matmul(transpose(c),matmul(SigTp,c)) Sigp = matmul(transpose(c),matmul(Sigp,c))
! Compute new density matrix in the AO basis ! Compute new density matrix in the AO basis
@ -311,7 +307,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Print results ! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigTp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) call print_qsGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
enddo enddo
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -334,7 +330,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis) deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,error,error_diis,F_diis)
! Perform BSE calculation ! Perform BSE calculation

3
src/GW/G0W0.f90
View File

@ -144,8 +144,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
else else
call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
end if end if

72
src/GW/GW_renormalization_factor.f90
View File

@ -1,72 +0,0 @@
subroutine GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute renormalization factor for GW
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
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) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: p,i,a,jb
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
! static COHSEX approximation
if(COHSEX) then
Z(:) = 1d0
return
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end if
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine

49
src/GW/GW_self_energy.f90
View File

@ -1,6 +1,6 @@
subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy and the renormalization factor
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -24,16 +24,18 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
integer :: i,j,a,b integer :: i,j,a,b
integer :: p,q,r integer :: p,q,r
integer :: jb integer :: jb
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(out) :: EcGM double precision,intent(out) :: EcGM
double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: Sig(nBas,nBas)
double precision,intent(out) :: Z(nBas)
! Initialize ! Initialize
SigC(:,:) = 0d0 Sig(:,:) = 0d0
Z(:) = 0d0
!-----------------------------! !-----------------------------!
! COHSEX static approximation ! ! COHSEX static approximation !
@ -47,7 +49,7 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do jb=1,nS do jb=1,nS
SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) Sig(p,q) = Sig(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb)
end do end do
end do end do
end do end do
@ -59,7 +61,7 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do r=nC+1,nBas-nR do r=nC+1,nBas-nR
do jb=1,nS do jb=1,nS
SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) Sig(p,q) = Sig(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb)
end do end do
end do end do
end do end do
@ -67,9 +69,11 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
EcGM = 0d0 EcGM = 0d0
do i=nC+1,nO do i=nC+1,nO
EcGM = EcGM + 0.5d0*SigC(i,i) EcGM = EcGM + 0.5d0*Sig(i,i)
end do end do
Z(:) = 0d0
else else
!----------------! !----------------!
@ -79,16 +83,20 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,i,q,p,eps) & !$OMP PRIVATE(jb,i,q,p,eps,num) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do jb=1,nS do jb=1,nS
do i=nC+1,nO do i=nC+1,nO
eps = e(p) - e(i) + Omega(jb) eps = e(p) - e(i) + Omega(jb)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2) num = 2d0*rho(p,i,jb)*rho(q,i,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -99,16 +107,20 @@ end do
! Virtual part of the correlation self-energy ! Virtual part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,a,q,p,eps) & !$OMP PRIVATE(jb,a,q,p,eps,num) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do jb=1,nS do jb=1,nS
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
eps = e(p) - e(a) - Omega(jb) eps = e(p) - e(a) - Omega(jb)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2) num = 2d0*rho(p,a,jb)*rho(q,a,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -122,12 +134,19 @@ end do
do jb=1,nS do jb=1,nS
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
eps = e(a) - e(i) + Omega(jb) eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) num = 4d0*rho(a,i,jb)*rho(a,i,jb)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
end do end do
end do end do
end if end if
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

39
src/GW/GW_self_energy_diag.f90
View File

@ -1,6 +1,6 @@
subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the self-energy ! Compute diagonal of the correlation part of the self-energy and the renormalization factor
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -22,16 +22,18 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
! Local variables ! Local variables
integer :: i,a,p,q,jb integer :: i,a,p,q,jb
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(out) :: SigC(nBas) double precision,intent(out) :: Sig(nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: EcGM double precision,intent(out) :: EcGM
! Initialize ! Initialize
SigC(:) = 0d0 Sig(:) = 0d0
Z(:) = 0d0
!----------------------------- !-----------------------------
! COHSEX static self-energy ! COHSEX static self-energy
@ -44,7 +46,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do jb=1,nS do jb=1,nS
SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/Omega(jb) Sig(p) = Sig(p) + 4d0*rho(p,i,jb)**2/Omega(jb)
end do end do
end do end do
end do end do
@ -54,7 +56,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do q=nC+1,nBas-nR do q=nC+1,nBas-nR
do jb=1,nS do jb=1,nS
SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb) Sig(p) = Sig(p) - 2d0*rho(p,q,jb)**2/Omega(jb)
end do end do
end do end do
end do end do
@ -63,7 +65,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
EcGM = 0d0 EcGM = 0d0
do i=nC+1,nO do i=nC+1,nO
EcGM = EcGM - SigC(i) EcGM = EcGM - Sig(i)
end do end do
!----------------------------- !-----------------------------
@ -77,8 +79,12 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do i=nC+1,nO do i=nC+1,nO
do jb=1,nS do jb=1,nS
eps = e(p) - e(i) + Omega(jb) eps = e(p) - e(i) + Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) num = 2d0*rho(p,i,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -88,8 +94,12 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do jb=1,nS do jb=1,nS
eps = e(p) - e(a) - Omega(jb) eps = e(p) - e(a) - Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) num = 2d0*rho(p,a,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
end do end do
@ -100,12 +110,19 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
do i=nC+1,nO do i=nC+1,nO
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do jb=1,nS do jb=1,nS
eps = e(a) - e(i) + Omega(jb) eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) num = 4d0*rho(a,i,jb)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
end do end do
end do end do
end if end if
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

3
src/GW/evGW.f90
View File

@ -167,8 +167,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
else else
call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
end if end if

3
src/GW/qsGW.f90
View File

@ -192,8 +192,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
else else
call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
endif endif