From 6cc9802d00e573ca85f1678cdca7e696c9afce5d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 21 Jul 2023 11:18:00 +0200 Subject: [PATCH] include B term in phBSE2@GW --- input/options | 4 +- src/GW/GW_excitation_density.f90 | 6 +- src/GW/GW_phBSE2_dynamic_kernel_A.f90 | 21 +++--- src/GW/GW_phBSE2_dynamic_kernel_B.f90 | 88 ++++++++++++++++++++++++ src/GW/GW_phBSE_dynamic_kernel_A.f90 | 7 +- src/GW/GW_phBSE_dynamic_kernel_B.f90 | 7 +- src/GW/GW_phBSE_dynamic_perturbation.f90 | 4 +- src/GW/GW_phBSE_static_kernel_A.f90 | 7 +- 8 files changed, 129 insertions(+), 15 deletions(-) create mode 100644 src/GW/GW_phBSE2_dynamic_kernel_B.f90 diff --git a/input/options b/input/options index 9db91ed..89427ee 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + F T F T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX TDA_W reg @@ -15,4 +15,4 @@ # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA - T F F T T + T T F T T diff --git a/src/GW/GW_excitation_density.f90 b/src/GW/GW_excitation_density.f90 index b00f3ca..a6ea714 100644 --- a/src/GW/GW_excitation_density.f90 +++ b/src/GW/GW_excitation_density.f90 @@ -6,7 +6,11 @@ subroutine GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) ! Input variables - integer,intent(in) :: nBas,nC,nO,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: XpY(nS,nS) diff --git a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 index 1d2e4d8..64cadc6 100644 --- a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn) +subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,KA_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -7,7 +7,12 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: W(nBas,nBas,nBas,nBas) double precision,intent(in) :: eGW(nBas) @@ -21,14 +26,14 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, ! Output variables - double precision,intent(out) :: A_dyn(nS,nS) + double precision,intent(out) :: KA_dyn(nS,nS) double precision,intent(out) :: ZA_dyn(nS,nS) ! Build dynamic A matrix jb = 0 -!$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmBSE,W,num,dem,eGW,nO,nBas,eta,nC,nR) +!$omp parallel do default(private) shared(KA_dyn,ZA_dyn,OmBSE,W,num,dem,eGW,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -44,13 +49,13 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, dem = OmBSE - eGW(a) - eGW(c) + eGW(k) + eGW(j) num = 2d0*W(j,k,i,c)*W(a,c,b,k) - A_dyn(ia,jb) = A_dyn(ia,jb) - num*dem/(dem**2 + eta**2) + KA_dyn(ia,jb) = KA_dyn(ia,jb) - num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 dem = OmBSE + eGW(i) - eGW(c) + eGW(k) - eGW(b) num = 2d0*W(j,c,i,k)*W(a,k,b,c) - A_dyn(ia,jb) = A_dyn(ia,jb) - num*dem/(dem**2 + eta**2) + KA_dyn(ia,jb) = KA_dyn(ia,jb) - num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 enddo @@ -62,7 +67,7 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, dem = OmBSE + eGW(i) + eGW(j) - eGW(c) - eGW(d) num = 2d0*W(a,j,c,d)*W(c,d,i,b) - A_dyn(ia,jb) = A_dyn(ia,jb) + num*dem/(dem**2 + eta**2) + KA_dyn(ia,jb) = KA_dyn(ia,jb) + num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 enddo @@ -74,7 +79,7 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, dem = OmBSE - eGW(a) - eGW(b) + eGW(k) + eGW(l) num = 2d0*W(a,j,k,l)*W(k,l,i,b) - A_dyn(ia,jb) = A_dyn(ia,jb) + num*dem/(dem**2 + eta**2) + KA_dyn(ia,jb) = KA_dyn(ia,jb) + num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do diff --git a/src/GW/GW_phBSE2_dynamic_kernel_B.f90 b/src/GW/GW_phBSE2_dynamic_kernel_B.f90 new file mode 100644 index 0000000..bc7d5b5 --- /dev/null +++ b/src/GW/GW_phBSE2_dynamic_kernel_B.f90 @@ -0,0 +1,88 @@ +subroutine GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,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 + double precision,intent(in) :: eta + double precision,intent(in) :: W(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGW(nBas) + +! Local variables + + double precision :: num,dem + double precision :: eps + integer :: i,j,a,b,ia,jb,kc,k,c,l,d + +! Output variables + + double precision,intent(out) :: KB_dyn(nS,nS) + +! Build dynamic B matrix + + jb = 0 + +!$omp parallel do default(private) shared(KB_dyn,W,num,dem,eGW,nO,nBas,eta,nC,nR) + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = (b-nO) + (j-1)*(nBas-nO) + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = (a-nO) + (i-1)*(nBas-nO) + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - eGW(a) - eGW(c) + eGW(k) + eGW(j) + num = 2d0*W(b,k,i,c)*W(a,c,j,k) + + KB_dyn(ia,jb) = KB_dyn(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + eGW(i) - eGW(c) + eGW(k) - eGW(b) + num = 2d0*W(b,c,i,k)*W(a,k,j,c) + + KB_dyn(ia,jb) = KB_dyn(ia,jb) - num*dem/(dem**2 + eta**2) + + enddo + enddo + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = + eGW(i) + eGW(j) - eGW(c) - eGW(d) + num = 2d0*W(a,b,c,d)*W(c,d,i,j) + + KB_dyn(ia,jb) = KB_dyn(ia,jb) + num*dem/(dem**2 + eta**2) + + enddo + enddo + + do k=nC+1,nO + do l=nC+1,nO + + dem = - eGW(a) - eGW(b) + eGW(k) + eGW(l) + num = 2d0*W(a,b,k,l)*W(k,l,i,j) + + KB_dyn(ia,jb) = KB_dyn(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + enddo + enddo + enddo + enddo +!$omp end parallel do + +end subroutine diff --git a/src/GW/GW_phBSE_dynamic_kernel_A.f90 b/src/GW/GW_phBSE_dynamic_kernel_A.f90 index 3020ab5..aaa528e 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90 @@ -7,7 +7,12 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) diff --git a/src/GW/GW_phBSE_dynamic_kernel_B.f90 b/src/GW/GW_phBSE_dynamic_kernel_B.f90 index 29365e3..afc5166 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_B.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_B.f90 @@ -7,7 +7,12 @@ subroutine GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) diff --git a/src/GW/GW_phBSE_dynamic_perturbation.f90 b/src/GW/GW_phBSE_dynamic_perturbation.f90 index fb2fc9b..09da1af 100644 --- a/src/GW/GW_phBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_phBSE_dynamic_perturbation.f90 @@ -65,7 +65,7 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e if(dophBSE2) then write(*,*) - write(*,*) '*** Second-order BSE static kernel activated! ***' + write(*,*) '*** Second-order BSE dynamic kernel activated! ***' write(*,*) allocate(W(nBas,nBas,nBas,nBas)) @@ -108,6 +108,8 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,-OmBSE(ia),KAm_dyn,ZAm_dyn) call GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,KB_dyn) + if(dophBSE2) call GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn) + ! Renormalization factor of the resonant and anti-resonant parts ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & diff --git a/src/GW/GW_phBSE_static_kernel_A.f90 b/src/GW/GW_phBSE_static_kernel_A.f90 index 2747173..fb3ec9d 100644 --- a/src/GW/GW_phBSE_static_kernel_A.f90 +++ b/src/GW/GW_phBSE_static_kernel_A.f90 @@ -7,7 +7,12 @@ subroutine GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KA ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)