From 2a668b4d4059310c263ce51159ee229f4bec4bcb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 4 Jul 2024 10:41:57 +0200 Subject: [PATCH] pp-dBSE beyond TDA --- src/GW/GW_phBSE_dynamic_perturbation.f90 | 28 ++++---- src/GW/GW_ppBSE_dynamic_kernel_B.f90 | 3 + src/GW/GW_ppBSE_dynamic_kernel_C.f90 | 3 + src/GW/GW_ppBSE_dynamic_kernel_D.f90 | 4 ++ src/GW/GW_ppBSE_dynamic_perturbation.f90 | 84 +++++++++++++++++------- 5 files changed, 84 insertions(+), 38 deletions(-) diff --git a/src/GW/GW_phBSE_dynamic_perturbation.f90 b/src/GW/GW_phBSE_dynamic_perturbation.f90 index 77ec532..06348be 100644 --- a/src/GW/GW_phBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_phBSE_dynamic_perturbation.f90 @@ -37,8 +37,8 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e integer :: maxS = 10 double precision :: gapGW - double precision,allocatable :: OmDyn(:) - double precision,allocatable :: ZDyn(:) + double precision,allocatable :: Om_dyn(:) + double precision,allocatable :: Z_dyn(:) double precision,allocatable :: X(:) 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 - 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)) 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 - ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) + Z_dyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) + Om_dyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) 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 - ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & - + dot_product(Y,matmul(ZAm_dyn,Y)) + Z_dyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) & - - dot_product(Y,matmul(KAm_dyn - KA_sta,Y)) & - + dot_product(X,matmul(KB_dyn - KB_sta,Y)) & - - dot_product(Y,matmul(KB_dyn - KB_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(X,matmul(KB_dyn - KB_sta,Y)) & + - dot_product(Y,matmul(KB_dyn - KB_sta,X)) end if - ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) - OmDyn(ia) = ZDyn(ia)*OmDyn(ia) + Z_dyn(ia) = 1d0/(1d0 - Z_dyn(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)') & - 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 write(*,*) '---------------------------------------------------------------------------------------------------' diff --git a/src/GW/GW_ppBSE_dynamic_kernel_B.f90 b/src/GW/GW_ppBSE_dynamic_kernel_B.f90 index d837bfe..8a7b4a4 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_B.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_B.f90 @@ -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) end do + KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) 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) end do + KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_kernel_C.f90 b/src/GW/GW_ppBSE_dynamic_kernel_C.f90 index c9f6445..e291f30 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_C.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_C.f90 @@ -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 end do + 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))) + 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) ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_kernel_D.f90 b/src/GW/GW_ppBSE_dynamic_kernel_D.f90 index a429a36..1497061 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_D.f90 @@ -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 end do + 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))) + 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 end do + KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl) ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_perturbation.f90 b/src/GW/GW_ppBSE_dynamic_perturbation.f90 index a72736d..aa7a68f 100644 --- a/src/GW/GW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_ppBSE_dynamic_perturbation.f90 @@ -38,15 +38,15 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO, ! Local variables - integer :: ab,ij + integer :: ab,ij,kl integer :: maxOO = 10 integer :: maxVV = 10 - double precision,allocatable :: Om1Dyn(:) - double precision,allocatable :: Om2Dyn(:) - double precision,allocatable :: Z1Dyn(:) - double precision,allocatable :: Z2Dyn(:) + double precision,allocatable :: Om1_dyn(:) + double precision,allocatable :: Om2_dyn(:) + double precision,allocatable :: Z1_dyn(:) + double precision,allocatable :: Z2_dyn(:) double precision,allocatable :: KB_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 - 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), & ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO)) - Om1Dyn(:) = 0d0 - Om2Dyn(:) = 0d0 - Z1Dyn(:) = 0d0 - Z2Dyn(:) = 0d0 + Om1_dyn(:) = 0d0 + Om2_dyn(:) = 0d0 + Z1_dyn(:) = 0d0 + Z2_dyn(:) = 0d0 if(dTDA) then 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) -! 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) - 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) + if(dTDA) then - Z1Dyn(ab) = dot_product(X1(ab,:),matmul(ZC_dyn,X1(ab,:))) - Om1Dyn(ab) = dot_product(X1(ab,:),matmul(KC_dyn - KC_sta,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) + + Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) + Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) - Z1Dyn(ab) = 1d0/(1d0 - Z1Dyn(ab)) - Om1Dyn(ab) = Z1Dyn(ab)*Om1Dyn(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)') & - 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 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(*,*) '---------------------------------------------------------------------------------------------------' + kl = 0 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) - 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) + if(dTDA) then - Z2Dyn(ij) = dot_product(Y2(ij,:),matmul(ZD_dyn,Y2(ij,:))) - Om2Dyn(ij) = dot_product(Y2(ij,:),matmul(KD_dyn - KD_sta,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) + + Z2_dyn(kl) = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) + Om2_dyn(kl) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) - Z2Dyn(ij) = 1d0/(1d0 - Z2Dyn(ij)) - Om2Dyn(ij) = Z2Dyn(ij)*Om2Dyn(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)') & - 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 write(*,*) '---------------------------------------------------------------------------------------------------'