From bb4729192ef9dbc6175826129fd9d990bfa26261 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 22 Aug 2024 15:23:08 +0200 Subject: [PATCH] GW_ppBSE_static_kernel_C for ispin = 1 --- src/GW/GW_ppBSE.f90 | 109 ++++++++++++++++++++++++---- src/GW/GW_ppBSE_static_kernel_C.f90 | 92 +++++++++++++++++------ src/GW/RG0W0.f90 | 56 +++++++++++++- 3 files changed, 221 insertions(+), 36 deletions(-) diff --git a/src/GW/GW_ppBSE.f90 b/src/GW/GW_ppBSE.f90 index c1d6daa..e26b2e1 100644 --- a/src/GW/GW_ppBSE.f90 +++ b/src/GW/GW_ppBSE.f90 @@ -66,6 +66,11 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, double precision,intent(out) :: EcBSE(nspin) + double precision :: t0, t1 + double precision :: tt0, tt1 + + call wall_time(t0) + !--------------------------------- ! Compute (singlet) RPA screening !--------------------------------- @@ -76,11 +81,25 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & Aph(nS,nS),Bph(nS,nS)) - call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) - if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call wall_time(tt0) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds' + call wall_time(tt0) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds' + + call wall_time(tt0) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds' + + call wall_time(tt0) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_excitation_density =',tt1-tt0,' seconds' deallocate(XpY_RPA,XmY_RPA,Aph,Bph) @@ -108,32 +127,63 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ! Compute BSE excitation energies - if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call wall_time(tt0) + call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_C =',tt1-tt0,' seconds' + call wall_time(tt0) + call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_D =',tt1-tt0,' seconds' + + call wall_time(tt0) + if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_B =',tt1-tt0,' seconds' + + call wall_time(tt0) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_C =',tt1-tt0,' seconds' + + call wall_time(tt0) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_D =',tt1-tt0,' seconds' + + call wall_time(tt0) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_B =',tt1-tt0,' seconds' Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + call wall_time(tt0) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR =',tt1-tt0,' seconds' + call wall_time(tt0) call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_transition_vectors =',tt1-tt0,' seconds' !----------------------------------------------------! ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! + call wall_time(tt0) if(dBSE) & call GW_ppBSE_dynamic_perturbation(ispin,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) - write(*,*) "Deallocate not done" + + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_dynamic_perturbation =',tt1-tt0,' seconds' + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) - write(*,*) "Deallocate done" end if !------------------- @@ -160,33 +210,66 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ! Compute BSE excitation energies + call wall_time(tt0) + call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_C =',tt1-tt0,' seconds' + + call wall_time(tt0) + call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_D =',tt1-tt0,' seconds' + + call wall_time(tt0) if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_B =',tt1-tt0,' seconds' + call wall_time(tt0) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_C =',tt1-tt0,' seconds' + call wall_time(tt0) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_D =',tt1-tt0,' seconds' + + call wall_time(tt0) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_B =',tt1-tt0,' seconds' Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + call wall_time(tt0) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR =',tt1-tt0,' seconds' + call wall_time(tt0) call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_transition_vectors =',tt1-tt0,' seconds' !----------------------------------------------------! ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! + call wall_time(tt0) if(dBSE) & call GW_ppBSE_dynamic_perturbation(ispin,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) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_dynamic_perturbation =',tt1-tt0,' seconds' deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) end if + call wall_time(t1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE =',t1-t0,' seconds' + end subroutine diff --git a/src/GW/GW_ppBSE_static_kernel_C.f90 b/src/GW/GW_ppBSE_static_kernel_C.f90 index ef21825..dfc9c75 100644 --- a/src/GW/GW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GW_ppBSE_static_kernel_C.f90 @@ -26,44 +26,94 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI double precision,external :: Kronecker_delta double precision :: chi double precision :: eps + double precision :: tmp_ab, lambda4, eta2 integer :: a,b,c,d,ab,cd,m + integer :: a0, aa + + double precision, allocatable :: Om_tmp(:) ! Output variables double precision,intent(out) :: KC(nVV,nVV) -! Initialization - - KC(:,:) = 0d0 - !---------------! ! Singlet block ! !---------------! if(ispin == 1) then - ab = 0 - do a=nO+1,nBas-nR - do b=a,nBas-nR - ab = ab + 1 + a0 = nBas - nR - nO + lambda4 = 4.d0 * lambda + eta2 = eta * eta + + allocate(Om_tmp(nS)) + + !$OMP PARALLEL DEFAULT(NONE) PRIVATE(m) SHARED(nS, eta2, Om, Om_tmp) + !$OMP DO + do m = 1, nS + Om_tmp(m) = Om(m) / (Om(m)*Om(m) + eta2) + enddo + !$OMP END DO + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, tmp_ab) & + !$OMP SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC) + !$OMP DO + do a = nO+1, nBas-nR + aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO + do b = a, nBas-nR + ab = aa + b + + tmp_ab = lambda4 + if(a .eq. b) then + tmp_ab = 0.7071067811865475d0 * lambda4 + endif + cd = 0 - do c=nO+1,nBas-nR - do d=c,nBas-nR + do c = nO + 1, nBas-nR + do d = c, nBas-nR cd = cd + 1 - chi = 0d0 - do m=1,nS - eps = Om(m)**2 + eta**2 - chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & - - rho(a,d,m)*rho(b,c,m)*Om(m)/eps - end do + KC(ab,cd) = 0d0 + do m = 1, nS + KC(ab,cd) = KC(ab,cd) - rho(a,c,m) * rho(b,d,m) * Om_tmp(m) & + - rho(a,d,m) * rho(b,c,m) * Om_tmp(m) + end do - KC(ab,cd) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC(ab,cd) = tmp_ab * KC(ab,cd) + if(c .eq. d) then + KC(ab,cd) = 0.7071067811865475d0 * KC(ab,cd) + endif + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL - end do - end do - end do - end do +! ab = 0 +! do a=nO+1,nBas-nR +! do b=a,nBas-nR +! ab = ab + 1 +! cd = 0 +! do c=nO+1,nBas-nR +! do d=c,nBas-nR +! cd = cd + 1 +! +! chi = 0d0 +! do m=1,nS +! eps = Om(m)**2 + eta**2 +! chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & +! - rho(a,d,m)*rho(b,c,m)*Om(m)/eps +! end do +! +! KC(ab,cd) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) +! +! end do +! end do +! end do +! end do end if diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index fc30b5d..1bd18e0 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -59,6 +59,11 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA double precision,allocatable :: eGWlin(:) double precision,allocatable :: eGW(:) + double precision :: t0, t1 + double precision :: tt0, tt1 + + call wall_time(t0) + ! Output variables ! Hello world @@ -101,26 +106,48 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute screening ! !-------------------! - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call wall_time(tt0) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds' + + call wall_time(tt0) if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds' + call wall_time(tt0) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds' + call wall_time(tt0) if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for print_excitation_energies =',tt1-tt0,' seconds' !--------------------------! ! Compute spectral weights ! !--------------------------! + call wall_time(tt0) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_excitation_density =',tt1-tt0,' seconds' !------------------------! ! Compute GW self-energy ! !------------------------! + call wall_time(tt0) if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_regularization =',tt1-tt0,' seconds' + call wall_time(tt0) call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_self_energy_diag =',tt1-tt0,' seconds' !-----------------------------------! ! Solve the quasi-particle equation ! @@ -128,6 +155,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Linearized or graphical solution? + call wall_time(tt0) eGWlin(:) = eHF(:) + Z(:)*SigC(:) if(linearize) then @@ -145,6 +173,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) end if + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for QP =',tt1-tt0,' seconds' ! Plot self-energy, renormalization factor, and spectral function @@ -158,19 +188,33 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute the RPA correlation energy + call wall_time(tt0) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds' + call wall_time(tt0) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds' + + call wall_time(tt0) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds' !--------------! ! Dump results ! !--------------! + call wall_time(tt0) call print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for print_RG0W0 =',tt1-tt0,' seconds' ! Perform BSE calculation + call wall_time(tt0) if(dophBSE) then call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) @@ -221,7 +265,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA end if end if + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for phBSE =',tt1-tt0,' seconds' + call wall_time(tt0) if(doppBSE) then call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) @@ -238,6 +285,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) end if + call wall_time(tt1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppBSE =',tt1-tt0,' seconds' ! end if @@ -251,4 +300,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA end if + call wall_time(t1) + write(*,'(A65,1X,F9.3,A8)') 'Wall time for RG0W0 =',t1-t0,' seconds' + end subroutine