From 2ecc80a2786e673be19fb1477a74e27ed2d8423a Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 24 Jul 2023 22:34:52 +0200 Subject: [PATCH] dynamic ppBSE@GW --- src/GW/G0W0.f90 | 2 +- src/GW/GW_phBSE.f90 | 6 +- src/GW/GW_ppBSE.f90 | 20 +++- src/GW/GW_ppBSE_dynamic_perturbation.f90 | 116 +++++++++++++++++++++++ src/GW/evGW.f90 | 2 +- src/GW/qsGW.f90 | 2 +- 6 files changed, 141 insertions(+), 7 deletions(-) create mode 100644 src/GW/GW_ppBSE_dynamic_perturbation.f90 diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 29a026f..c630726 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -226,7 +226,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90 index 8da8502..b74828e 100644 --- a/src/GW/GW_phBSE.f90 +++ b/src/GW/GW_phBSE.f90 @@ -120,9 +120,9 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO, 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) - !------------------------------------------------- - ! Compute the dynamical screening at the BSE level - !------------------------------------------------- + !----------------------------------------------------! + ! Compute the dynamical screening at the phBSE level ! + !----------------------------------------------------! if(dBSE) & call GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & diff --git a/src/GW/GW_ppBSE.f90 b/src/GW/GW_ppBSE.f90 index 7ab6990..a8497b1 100644 --- a/src/GW/GW_ppBSE.f90 +++ b/src/GW/GW_ppBSE.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies at the pp level @@ -9,6 +9,8 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -122,6 +124,14 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole call print_transition_vectors_pp(.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 GW_ppBSE_dynamic_perturbation(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 if @@ -167,6 +177,14 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole call print_transition_vectors_pp(.false.,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 GW_ppBSE_dynamic_perturbation(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 if diff --git a/src/GW/GW_ppBSE_dynamic_perturbation.f90 b/src/GW/GW_ppBSE_dynamic_perturbation.f90 new file mode 100644 index 0000000..1634f8d --- /dev/null +++ b/src/GW/GW_ppBSE_dynamic_perturbation.f90 @@ -0,0 +1,116 @@ +subroutine GW_ppBSE_dynamic_perturbation(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) + +! Compute dynamical effects via perturbation theory for BSE + + implicit none + include 'parameters.h' + +! Input variables + + 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 + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + double precision,intent(in) :: KB_sta(nVV,nOO) + double precision,intent(in) :: KC_sta(nVV,nVV) + double precision,intent(in) :: KD_sta(nOO,nOO) + +! Local variables + + integer :: ab,ij + + 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 :: KB_dyn(:,:) + double precision,allocatable :: KC_dyn(:,:) + double precision,allocatable :: KD_dyn(:,:) + double precision,allocatable :: ZC_dyn(:,:) + double precision,allocatable :: ZD_dyn(:,:) + +! Memory allocation + + allocate(Om1Dyn(nVV),Om2Dyn(nOO),Z1Dyn(maxVV),Z2Dyn(maxOO), & + KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), & + ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO)) + + if(dTDA) then + write(*,*) + write(*,*) '*** dynamical TDA activated ***' + write(*,*) + end if + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static ppBSE double electron attachment energies ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + 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(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))) + Om1Dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) + + Z1Dyn(ab) = 1d0/(1d0 - Z1Dyn(ab)) + Om1Dyn(ab) = Z1Dyn(ab)*Om1Dyn(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) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static ppBSE double electron detachment energies ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do ij=1,min(nOO,maxOO) + +! 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(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))) + Om2Dyn(ij) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) + + Z2Dyn(ij) = 1d0/(1d0 - Z2Dyn(ij)) + Om2Dyn(ij) = Z2Dyn(ij)*Om2Dyn(ij) + + 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) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index a763a90..e333956 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -276,7 +276,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 50bbe16..e031da6 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -341,7 +341,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) + call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------'