LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-31 08:36:05 +01:00
Code Releases Activity

ppBSE@GW@GHF

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-06 13:47:15 +02:00
parent 3111666d15
commit fa69888f48
3 changed files with 144 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
mol/LiF.xyz
View File

@ -1,4 +1,4 @@
2
Li 0.0000 0.0000 0.0000
F 0.0000 0.0000 1.5783
F 0.0000 0.0000 1.58753

33
src/GW/GG0W0.f90
View File

@ -1,5 +1,5 @@
subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Perform G0W0 calculation
implicit none
@ -31,7 +31,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: EGHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas)
@ -156,7 +156,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Dump results !
!--------------!
call print_GG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
call print_GG0W0(nBas,nO,eHF,ENuc,EGHF,SigC,Z,eGW,EcRPA,EcGM)
! Deallocate memory
@ -171,7 +171,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF correlation energy = ',EcBSE,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy = ',ENuc + ERHF + EcBSE,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -198,7 +198,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + EGHF + EcBSE(1) + EcBSE(2),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
@ -206,26 +206,25 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
end if
! if(doppBSE) then
if(doppBSE) then
! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
call GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 correlation energy =',EcBSE,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 total energy =',ENuc + EGHF + EcBSE,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! end if
end if
! Testing zone
if(dotest) then
call dump_test_value('G','G0W0 correlation energy',EcRPA)
call dump_test_value('G','RPA@G0W0 correlation energy',EcRPA)
call dump_test_value('G','Tr@ppBSE@G0W0 correlation energy',EcBSE)
call dump_test_value('G','G0W0 HOMO energy',eGW(nO))
call dump_test_value('G','G0W0 LUMO energy',eGW(nO+1))

127
src/GW/GGW_ppBSE.f90 Normal file
View File

@ -0,0 +1,127 @@
subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
! Compute the Bethe-Salpeter excitation energies at the pp level based on a GHF reference
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
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) :: eW(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: isp_W
logical :: dRPA = .false.
logical :: dRPA_W = .true.
integer :: nOO
integer :: nVV
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Cpp(:,:)
double precision,allocatable :: Dpp(:,:)
double precision,allocatable :: Om1(:)
double precision,allocatable :: X1(:,:)
double precision,allocatable :: Y1(:,:)
double precision,allocatable :: Om2(:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: Y2(:,:)
double precision,allocatable :: KB_sta(:,:)
double precision,allocatable :: KC_sta(:,:)
double precision,allocatable :: KD_sta(:,:)
! Output variables
double precision,intent(out) :: EcBSE
!-----------------------!
! Compute RPA screening !
!-----------------------!
isp_W = 3
EcRPA = 0d0
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
Aph(nS,nS),Bph(nS,nS))
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
ispin = 4
EcBSE = 0d0
nOO = nO*(nO-1)/2
nVV = nV*(nV-1)/2
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
! Compute BSE excitation energies
call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
!----------------------------------------------------!
! Compute the dynamical screening at the ppBSE level !
!----------------------------------------------------!
! if(dBSE) &
! call GGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
end subroutine