PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

pp-dBSE beyond TDA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-07-04 10:41:57 +02:00
parent 38a83e4a1d
commit 2a668b4d40
5 changed files with 84 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/GW/GW_phBSE_dynamic_perturbation.f90
View File

@ -37,8 +37,8 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
integer :: maxS = 10 integer :: maxS = 10
double precision :: gapGW double precision :: gapGW
double precision,allocatable :: OmDyn(:) double precision,allocatable :: Om_dyn(:)
double precision,allocatable :: ZDyn(:) double precision,allocatable :: Z_dyn(:)
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
@ -53,7 +53,7 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
! Memory allocation ! Memory allocation
allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS), & allocate(Om_dyn(maxS),Z_dyn(maxS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS), &
KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS)) KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS))
if(dTDA) then if(dTDA) then
@ -96,8 +96,8 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
if(dTDA) then if(dTDA) then
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) Z_dyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) Om_dyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X))
else else
@ -110,21 +110,21 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
! Renormalization factor of the resonant and anti-resonant parts ! Renormalization factor of the resonant and anti-resonant parts
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & Z_dyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &
+ dot_product(Y,matmul(ZAm_dyn,Y)) + dot_product(Y,matmul(ZAm_dyn,Y))
OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) & Om_dyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) &
- dot_product(Y,matmul(KAm_dyn - KA_sta,Y)) & - dot_product(Y,matmul(KAm_dyn - KA_sta,Y)) &
+ dot_product(X,matmul(KB_dyn - KB_sta,Y)) & + dot_product(X,matmul(KB_dyn - KB_sta,Y)) &
- dot_product(Y,matmul(KB_dyn - KB_sta,X)) - dot_product(Y,matmul(KB_dyn - KB_sta,X))
end if end if
ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) Z_dyn(ia) = 1d0/(1d0 - Z_dyn(ia))
OmDyn(ia) = ZDyn(ia)*OmDyn(ia) Om_dyn(ia) = Z_dyn(ia)*Om_dyn(ia)
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+Om_dyn(ia))*HaToeV,Om_dyn(ia)*HaToeV,Z_dyn(ia)
end do end do
write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'

3
src/GW/GW_ppBSE_dynamic_kernel_B.f90
View File

@ -76,6 +76,7 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
end do end do
KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
end do end do
@ -121,7 +122,9 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
end do end do
KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij) KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)
end do end do
end do end do

3
src/GW/GW_ppBSE_dynamic_kernel_C.f90
View File

@ -80,8 +80,10 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do end do
KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
end do end do
end do end do
@ -132,6 +134,7 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG
KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd) KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)
ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd) ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)
end do end do
end do end do

4
src/GW/GW_ppBSE_dynamic_kernel_D.f90
View File

@ -81,8 +81,10 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do end do
KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
end do end do
end do end do
@ -130,8 +132,10 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do end do
KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl) KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl)
ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl) ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl)
end do end do
end do end do

84
src/GW/GW_ppBSE_dynamic_perturbation.f90
View File

@ -38,15 +38,15 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
! Local variables ! Local variables
integer :: ab,ij integer :: ab,ij,kl
integer :: maxOO = 10 integer :: maxOO = 10
integer :: maxVV = 10 integer :: maxVV = 10
double precision,allocatable :: Om1Dyn(:) double precision,allocatable :: Om1_dyn(:)
double precision,allocatable :: Om2Dyn(:) double precision,allocatable :: Om2_dyn(:)
double precision,allocatable :: Z1Dyn(:) double precision,allocatable :: Z1_dyn(:)
double precision,allocatable :: Z2Dyn(:) double precision,allocatable :: Z2_dyn(:)
double precision,allocatable :: KB_dyn(:,:) double precision,allocatable :: KB_dyn(:,:)
double precision,allocatable :: KC_dyn(:,:) double precision,allocatable :: KC_dyn(:,:)
@ -56,14 +56,14 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
! Memory allocation ! Memory allocation
allocate(Om1Dyn(nOO),Om2Dyn(nOO),Z1Dyn(nOO),Z2Dyn(nOO), & allocate(Om1_dyn(nOO),Om2_dyn(nOO),Z1_dyn(nOO),Z2_dyn(nOO), &
KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), & KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), &
ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO)) ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO))
Om1Dyn(:) = 0d0 Om1_dyn(:) = 0d0
Om2Dyn(:) = 0d0 Om2_dyn(:) = 0d0
Z1Dyn(:) = 0d0 Z1_dyn(:) = 0d0
Z2Dyn(:) = 0d0 Z2_dyn(:) = 0d0
if(dTDA) then if(dTDA) then
write(*,*) write(*,*)
@ -79,17 +79,34 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
do ab=1,min(nVV,maxVV) do ab=1,min(nVV,maxVV)
! if(.not.dTDA) call GW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,OmBSE(ab),KB_dyn) if(dTDA) then
call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
Z1Dyn(ab) = dot_product(X1(ab,:),matmul(ZC_dyn,X1(ab,:))) call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
Om1Dyn(ab) = dot_product(X1(ab,:),matmul(KC_dyn - KC_sta,X1(ab,:)))
Z1Dyn(ab) = 1d0/(1d0 - Z1Dyn(ab)) Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
Om1Dyn(ab) = Z1Dyn(ab)*Om1Dyn(ab) Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
else
call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KB_dyn)
call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KD_dyn,ZD_dyn)
Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) &
+ dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab)))
Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) &
- dot_product(Y1(:,ab),matmul(KD_dyn - KD_sta,Y1(:,ab))) &
+ dot_product(X1(:,ab),matmul(KB_dyn - KB_sta,Y1(:,ab))) &
- dot_product(Y1(:,ab),matmul(transpose(KB_dyn - KB_sta),X1(:,ab)))
end if
Z1_dyn(ab) = 1d0/(1d0 - Z1_dyn(ab))
Om1_dyn(ab) = Z1_dyn(ab)*Om1_dyn(ab)
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ab,Om1(ab)*HaToeV,(Om1(ab)+Om1Dyn(ab))*HaToeV,Om1Dyn(ab)*HaToeV,Z1Dyn(ab) ab,Om1(ab)*HaToeV,(Om1(ab)+Om1_dyn(ab))*HaToeV,Om1_dyn(ab)*HaToeV,Z1_dyn(ab)
end do end do
write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'
@ -101,19 +118,38 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'
kl = 0
do ij=max(1,nOO+1-maxOO),nOO do ij=max(1,nOO+1-maxOO),nOO
kl = kl + 1
! if(.not.dTDA) call GW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,OmBSE(ab),KB_dyn) if(dTDA) then
call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn)
Z2Dyn(ij) = dot_product(Y2(ij,:),matmul(ZD_dyn,Y2(ij,:))) call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KD_dyn,ZD_dyn)
Om2Dyn(ij) = dot_product(Y2(ij,:),matmul(KD_dyn - KD_sta,Y2(ij,:)))
Z2Dyn(ij) = 1d0/(1d0 - Z2Dyn(ij)) Z2_dyn(kl) = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
Om2Dyn(ij) = Z2Dyn(ij)*Om2Dyn(ij) Om2_dyn(kl) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
else
call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KB_dyn)
call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KC_dyn,ZC_dyn)
call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KD_dyn,ZD_dyn)
Z2_dyn(kl) = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) &
+ dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
Om2_dyn(kl) = dot_product(X2(:,ij),matmul(KC_dyn - KC_sta,X2(:,ij))) &
- dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) &
+ dot_product(X2(:,ij),matmul(KB_dyn - KB_sta,Y2(:,ij))) &
- dot_product(Y2(:,ij),matmul(transpose(KB_dyn - KB_sta),X2(:,ij)))
end if
Z2_dyn(kl) = 1d0/(1d0 - Z2_dyn(kl))
Om2_dyn(kl) = Z2_dyn(kl)*Om2_dyn(kl)
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ij,Om2(ij)*HaToeV,(Om2(ij)+Om2Dyn(ij))*HaToeV,Om2Dyn(ij)*HaToeV,Z2Dyn(ij) ij,Om2(ij)*HaToeV,(Om2(ij)+Om2_dyn(kl))*HaToeV,Om2_dyn(kl)*HaToeV,Z2_dyn(kl)
end do end do
write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'