From 61288e474d32e88907afecc99eedeebc60c8e017 Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 25 Oct 2024 14:44:35 +0200 Subject: [PATCH] ppBSE dynamic correction for the block D is working --- src/GW/GGW_ppBSE.f90 | 5 +-- src/GW/GGW_ppBSE_dynamic_kernel_D.f90 | 40 ++++++++++------------- src/GW/GGW_ppBSE_dynamic_perturbation.f90 | 4 +-- src/GW/GGW_ppBSE_static_kernel_D.f90 | 2 +- src/GW/RGW_ppBSE.f90 | 5 +-- src/GW/RGW_ppBSE_dynamic_kernel_D.f90 | 9 +++-- 6 files changed, 31 insertions(+), 34 deletions(-) diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90 index 5c48782..9303dff 100644 --- a/src/GW/GGW_ppBSE.f90 +++ b/src/GW/GGW_ppBSE.f90 @@ -72,7 +72,8 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int, if(.not.TDA_W) call phGLR_B(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phGLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - +! call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmRPA,XpY_RPA,XmY_RPA) + call GGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) deallocate(XpY_RPA,XmY_RPA,Aph,Bph) @@ -115,7 +116,7 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int, ! Upfolded ppBSE ! !----------------! - call GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) +! call GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) 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_D.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 index ba1a670..07d5e67 100644 --- a/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 @@ -52,33 +52,29 @@ subroutine GGW_ppBSE_dynamic_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om, kl = kl + 1 do m=1,nS + + num = (rho(i,k,m)*rho(j,l,m) - rho(j,k,m)*rho(i,l,m))/2 + dem = - OmBSE - Om(m) + eGW(j) + eGW(l) + 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(i,k,m)*rho(j,l,m) + dem = - OmBSE - Om(m) + eGW(i) + eGW(k) + 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 - 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(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 - - dem = OmBSE - eGW(l) + 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(l) + 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 - Om(m) + eGW(i) + eGW(l) + 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 - Om(m) + eGW(j) + eGW(k) + 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 + KD_dyn(ij,kl) = 0.5d0*KD_dyn(ij,kl) + ZD_dyn(ij,kl) = 0.5d0*ZD_dyn(ij,kl) + end do end do diff --git a/src/GW/GGW_ppBSE_dynamic_perturbation.f90 b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 index 344c207..968fab3 100644 --- a/src/GW/GGW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/GGW_ppBSE_dynamic_perturbation.f90 @@ -39,8 +39,8 @@ subroutine GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,e integer :: ab,ij,kl - integer :: maxOO = 0 - integer :: maxVV = 10 + integer :: maxOO = 20 + integer :: maxVV = 0 double precision,allocatable :: Om1_dyn(:) double precision,allocatable :: Om2_dyn(:) diff --git a/src/GW/GGW_ppBSE_static_kernel_D.f90 b/src/GW/GGW_ppBSE_static_kernel_D.f90 index 57e5f10..37c0672 100644 --- a/src/GW/GGW_ppBSE_static_kernel_D.f90 +++ b/src/GW/GGW_ppBSE_static_kernel_D.f90 @@ -55,7 +55,7 @@ subroutine GGW_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,r + rho(i,l,m)*rho(j,k,m)*Om(m)/eps end do - KD(ij,kl) = lambda*chi + KD(ij,kl) = 2d0*lambda*chi end do end do diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index ba1c10a..baf1d8a 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -103,6 +103,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) +! call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmRPA,XpY_RPA,XmY_RPA) call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) @@ -257,7 +258,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS ! Upfolded ppBSE ! !----------------! - call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) + ! call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) deallocate(KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2) @@ -375,7 +376,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS ! Upfolded ppBSE ! !----------------! - call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) + ! call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) deallocate(KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2) diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 index bd06642..6e5e332 100644 --- a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 @@ -76,8 +76,8 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e end do - KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) - ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + KD_dyn(ij,kl) = KD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -101,7 +101,6 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e do m=1,nS num = (rho(i,k,m)*rho(j,l,m) - rho(j,k,m)*rho(i,l,m))/2 -! dem = - Om(m) dem = - OmBSE - Om(m) + eGW(j) + eGW(l) 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 @@ -120,8 +119,8 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e end do - KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl) - ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl) + KD_dyn(ij,kl) = KD_dyn(ij,kl) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) end do end do