LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-19 04:22:39 +01:00
Code Releases Activity

fix bug in ppBSE

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-21 12:00:08 +02:00
parent d11bebafe1
commit 1d9725546d
3 changed files with 46 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/options
View File

@ -15,4 +15,4 @@
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA
T T F T F F F T F F

32
src/GW/GW_phBSE.f90
View File

@ -45,8 +45,8 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
double precision,allocatable :: XpY_BSE(:,:) double precision,allocatable :: XpY_BSE(:,:)
double precision,allocatable :: XmY_BSE(:,:) double precision,allocatable :: XmY_BSE(:,:)
double precision,allocatable :: A_sta(:,:) double precision,allocatable :: Aph(:,:)
double precision,allocatable :: B_sta(:,:) double precision,allocatable :: Bph(:,:)
double precision,allocatable :: KA_sta(:,:) double precision,allocatable :: KA_sta(:,:)
double precision,allocatable :: KB_sta(:,:) double precision,allocatable :: KB_sta(:,:)
@ -60,7 +60,7 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
! Memory allocation ! Memory allocation
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
A_sta(nS,nS),B_sta(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), & Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), &
OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
!--------------------------------- !---------------------------------
@ -70,10 +70,10 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
isp_W = 1 isp_W = 1
EcRPA = 0d0 EcRPA = 0d0
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,A_sta) 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,B_sta) 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,A_sta,B_sta,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
@ -90,8 +90,8 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
! Compute BSE excitation energies ! Compute BSE excitation energies
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
! Second-order BSE static kernel ! Second-order BSE static kernel
@ -112,10 +112,10 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
end if end if
A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) Aph(:,:) = Aph(:,:) + KA_sta(:,:)
if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation('phBSE@GW ',ispin,nS,OmBSE) call print_excitation('phBSE@GW ',ispin,nS,OmBSE)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
@ -141,13 +141,13 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
! Compute BSE excitation energies ! Compute BSE excitation energies
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) Aph(:,:) = Aph(:,:) + KA_sta(:,:)
if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation('phBSE@GW ',ispin,nS,OmBSE) call print_excitation('phBSE@GW ',ispin,nS,OmBSE)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)

58
src/GW/GW_ppBSE.f90
View File

@ -35,11 +35,8 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
integer :: nOO integer :: nOO
integer :: nVV integer :: nVV
double precision,allocatable :: A_sta(:,:) double precision,allocatable :: Aph(:,:)
double precision,allocatable :: B_sta(:,:) double precision,allocatable :: Bph(:,:)
double precision,allocatable :: KA_sta(:,:)
double precision,allocatable :: KB_sta(:,:)
double precision :: EcRPA double precision :: EcRPA
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
@ -59,9 +56,9 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
double precision,allocatable :: X2(:,:) double precision,allocatable :: X2(:,:)
double precision,allocatable :: Y2(:,:) double precision,allocatable :: Y2(:,:)
double precision,allocatable :: WB(:,:) double precision,allocatable :: KB_sta(:,:)
double precision,allocatable :: WC(:,:) double precision,allocatable :: KC_sta(:,:)
double precision,allocatable :: WD(:,:) double precision,allocatable :: KD_sta(:,:)
! Output variables ! Output variables
@ -74,14 +71,16 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
isp_W = 1 isp_W = 1
EcRPA = 0d0 EcRPA = 0d0
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,A_sta) allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) Aph(nS,nS),Bph(nS,nS))
call phLR(TDA_W,nS,A_sta,B_sta,EcRPA,OmRPA,XpY_RPA,XmY_RPA) 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 GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
!------------------- !-------------------
! Singlet manifold ! Singlet manifold
@ -95,6 +94,7 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
write(*,*) write(*,*)
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0
nOO = nO*(nO+1)/2 nOO = nO*(nO+1)/2
nVV = nV*(nV+1)/2 nVV = nV*(nV+1)/2
@ -102,27 +102,27 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
! Compute BSE excitation energies ! Compute BSE excitation energies
if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC) call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD) call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KD_sta)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) 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) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) + WB(:,:) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + WC(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + WD(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,WB,WC,WD) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
end if end if
@ -146,28 +146,28 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
! Compute BSE excitation energies ! Compute BSE excitation energies
if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC) call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD) call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KD_sta)
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) 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) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) + WB(:,:) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + WC(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + WD(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,WB,WC,WD) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
end if end if