From d8678bfe4f3900c046c2d42cd75c4dd53d81a2f7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 21 Jul 2023 10:45:10 +0200 Subject: [PATCH] fix bug in dBSE@GW --- input/options | 4 ++-- src/GW/G0W0.f90 | 1 + src/GW/GW_phBSE.f90 | 10 +++++----- src/GW/GW_phBSE_dynamic_kernel_A.f90 | 20 +++++--------------- src/QuAcK/QuAcK.f90 | 2 +- 5 files changed, 14 insertions(+), 23 deletions(-) diff --git a/input/options b/input/options index dfa37b1..773395e 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 T T T F + T F F T T diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 8a12163..da07811 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -183,6 +183,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) + if(exchange_kernel) then EcBSE(1) = 0.5d0*EcBSE(1) diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90 index afc48c6..ad3b0a9 100644 --- a/src/GW/GW_phBSE.f90 +++ b/src/GW/GW_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine GW_phBSE(dophBSE2,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 @@ -7,7 +7,7 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,n ! Input variables - logical,intent(in) :: BSE2 + logical,intent(in) :: dophBSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -95,7 +95,7 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,n ! Second-order BSE static kernel - if(BSE2) then + if(dophBSE2) then allocate(W(nBas,nBas,nBas,nBas)) @@ -125,7 +125,7 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,n !------------------------------------------------- if(dBSE) & - call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) end if @@ -157,7 +157,7 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,n !------------------------------------------------- if(dBSE) & - call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) end if diff --git a/src/GW/GW_phBSE_dynamic_kernel_A.f90 b/src/GW/GW_phBSE_dynamic_kernel_A.f90 index ea5a244..3020ab5 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn) +subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,KA_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -23,18 +23,18 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh ! 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) ! Initialization - A_dyn(:,:) = 0d0 + KA_dyn(:,:) = 0d0 ZA_dyn(:,:) = 0d0 ! Build dynamic A matrix jb = 0 -!$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmRPA,OmBSE,eGW,rho_RPA,nO,nBas,nS,chi,eps,eta,nC,nR,lambda) +!$omp parallel do default(private) shared(KA_dyn,ZA_dyn,OmRPA,OmBSE,eGW,rho_RPA,nO,nBas,nS,chi,eps,eta,nC,nR,lambda) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -45,16 +45,6 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh ia = (a-nO) + (i-1)*(nBas-nO) chi = 0d0 - - do kc=1,nS - eps = OmRPA(kc) - chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) - enddo - - A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi - - chi = 0d0 - do kc=1,nS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) @@ -65,7 +55,7 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh enddo - A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi + KA_dyn(ia,jb) = KA_dyn(ia,jb) - 2d0*lambda*chi chi = 0d0 do kc=1,nS diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 025b590..3a9f8e3 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -145,7 +145,7 @@ program QuAcK maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW,COHSEX,TDA_W, & maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dBSE,dTDA,doppBSE,dophBSE2) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA) !------------------------------------------------------------------------ ! Read input information