From 88b0424cd36b109193501dd98fd406342185f0d1 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 26 Sep 2024 14:26:52 +0200 Subject: [PATCH] GGW ppBSE dyn --- src/GW/GGW_ppBSE.f90 | 6 +- src/GW/GGW_ppBSE_dynamic_kernel_B.f90 | 82 ++++++++++++ src/GW/GGW_ppBSE_dynamic_kernel_C.f90 | 88 +++++++++++++ src/GW/GGW_ppBSE_dynamic_kernel_D.f90 | 88 +++++++++++++ src/GW/GGW_ppBSE_dynamic_perturbation.f90 | 152 ++++++++++++++++++++++ 5 files changed, 413 insertions(+), 3 deletions(-) create mode 100644 src/GW/GGW_ppBSE_dynamic_kernel_B.f90 create mode 100644 src/GW/GGW_ppBSE_dynamic_kernel_C.f90 create mode 100644 src/GW/GGW_ppBSE_dynamic_kernel_D.f90 create mode 100644 src/GW/GGW_ppBSE_dynamic_perturbation.f90 diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90 index 5adeba5..49ac42d 100644 --- a/src/GW/GGW_ppBSE.f90 +++ b/src/GW/GGW_ppBSE.f90 @@ -111,9 +111,9 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int, ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! -! if(dBSE) & -! call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,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) + if(dBSE) & + call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,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) diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 new file mode 100644 index 0000000..52dc5b4 --- /dev/null +++ b/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 @@ -0,0 +1,82 @@ +subroutine GGW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW,Om,rho,KB_dyn) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision,external :: Kronecker_delta + double precision :: dem,num + integer :: m + integer :: a,b,i,j + integer :: ab,ij + +! Output variables + + double precision,intent(out) :: KB_dyn(nVV,nOO) + +! Initialization + + KB_dyn(:,:) = 0d0 + +! Build dynamic B matrix + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + do m=1,nS + + dem = eGW(j) - Om(m) - eGW(b) + num = rho(a,i,m)*rho(b,j,m) + + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + + dem = eGW(j) - Om(m) - eGW(a) + num = rho(b,i,m)*rho(a,j,m) + + KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + + dem = eGW(i) - Om(m) - eGW(a) + num = rho(a,i,m)*rho(b,j,m) + + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + + dem = eGW(i) - Om(m) - eGW(b) + num = rho(b,i,m)*rho(a,j,m) + + KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + + end do + + end do + end do + + end do + end do + +end subroutine diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 new file mode 100644 index 0000000..76f02ce --- /dev/null +++ b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 @@ -0,0 +1,88 @@ +subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om,rho,OmBSE,KC_dyn,ZC_dyn) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision,external :: Kronecker_delta + double precision :: dem,num + integer :: m + integer :: a,b,c,d + integer :: ab,cd + +! Output variables + + double precision,intent(out) :: KC_dyn(nVV,nVV) + double precision,intent(out) :: ZC_dyn(nVV,nVV) + +! Initialization + + KC_dyn(:,:) = 0d0 + ZC_dyn(:,:) = 0d0 + +! Build dynamic C matrix + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + + do m=1,nS + + dem = OmBSE - eGW(c) - Om(m) - eGW(b) + num = rho(a,c,m)*rho(b,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(c) - Om(m) - eGW(a) + num = rho(b,c,m)*rho(a,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(d) - Om(m) - eGW(a) + num = rho(a,c,m)*rho(b,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(d) - Om(m) - eGW(b) + num = rho(b,c,m)*rho(a,d,m) + + KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + end do + + end do + end do + + end do + end do + +end subroutine diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 new file mode 100644 index 0000000..73dbe43 --- /dev/null +++ b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 @@ -0,0 +1,88 @@ +subroutine GGW_ppBSE_dynamic_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om,rho,OmBSE,KD_dyn,ZD_dyn) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + 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 + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision,external :: Kronecker_delta + double precision :: dem,num + integer :: m + integer :: i,j,k,l + integer :: ij,kl + +! Output variables + + double precision,intent(out) :: KD_dyn(nOO,nOO) + double precision,intent(out) :: ZD_dyn(nOO,nOO) + +! Initialization + + KD_dyn(:,:) = 0d0 + ZD_dyn(:,:) = 0d0 + +! Build dynamic D matrix + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + do m=1,nS + + dem = OmBSE - eGW(k) + Om(m) - eGW(i) + num = rho(i,k,m)*rho(j,l,m) + + KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(k) + Om(m) - eGW(j) + num = rho(j,k,m)*rho(i,l,m) + + KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(l) + Om(m) - eGW(j) + num = rho(i,k,m)*rho(j,l,m) + + KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - eGW(l) + Om(m) - eGW(i) + num = rho(j,k,m)*rho(i,l,m) + + KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + end do + + end do + end do + + end do + end do + +end subroutine diff --git a/src/GW/GGW_ppBSE_dynamic_perturbation.f90 b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 new file mode 100644 index 0000000..5e5f037 --- /dev/null +++ b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 @@ -0,0 +1,152 @@ +subroutine GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,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) :: nOrb + 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(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eW(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nOrb,nOrb,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,kl + + integer :: maxOO = 10 + integer :: maxVV = 0 + + 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(:,:) + double precision,allocatable :: KD_dyn(:,:) + double precision,allocatable :: ZC_dyn(:,:) + double precision,allocatable :: ZD_dyn(:,:) + +! Memory allocation + + allocate(Om1_dyn(maxVV),Om2_dyn(maxOO),Z1_dyn(maxVV),Z2_dyn(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(dTDA) then + + call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,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))) + + else + + call GGW_ppBSE_dynamic_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn) + call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn) + call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,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)+Om1_dyn(ab))*HaToeV,Om1_dyn(ab)*HaToeV,Z1_dyn(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(*,*) '---------------------------------------------------------------------------------------------------' + + kl = 0 + do ij=nOO,max(1,nOO+1-maxOO),-1 + kl = kl + 1 + + if(dTDA) then + + call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,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))) + + else + + call GGW_ppBSE_dynamic_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn) + call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KC_dyn,ZC_dyn) + call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,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)') & + kl,Om2(ij)*HaToeV,(Om2(ij)+Om2_dyn(kl))*HaToeV,Om2_dyn(kl)*HaToeV,Z2_dyn(kl) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine