4
1
mirror of https://github.com/pfloos/quack synced 2024-10-28 09:48:27 +01:00

ppBSE dynamic correction for the block D is working

This commit is contained in:
Antoine Marie 2024-10-25 14:44:35 +02:00
parent 7c972806a0
commit 61288e474d
6 changed files with 31 additions and 34 deletions

View File

@ -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) 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 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) call GGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
deallocate(XpY_RPA,XmY_RPA,Aph,Bph) 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 ! ! 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) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)

View File

@ -52,33 +52,29 @@ subroutine GGW_ppBSE_dynamic_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om,
kl = kl + 1 kl = kl + 1
do m=1,nS 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) dem = - OmBSE - Om(m) + eGW(i) + eGW(k)
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
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) dem = - OmBSE - Om(m) + eGW(i) + eGW(l)
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) dem = - OmBSE - Om(m) + eGW(j) + eGW(k)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**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
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
end do 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
end do end do

View File

@ -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 :: ab,ij,kl
integer :: maxOO = 0 integer :: maxOO = 20
integer :: maxVV = 10 integer :: maxVV = 0
double precision,allocatable :: Om1_dyn(:) double precision,allocatable :: Om1_dyn(:)
double precision,allocatable :: Om2_dyn(:) double precision,allocatable :: Om2_dyn(:)

View File

@ -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 + rho(i,l,m)*rho(j,k,m)*Om(m)/eps
end do end do
KD(ij,kl) = lambda*chi KD(ij,kl) = 2d0*lambda*chi
end do end do
end do end do

View File

@ -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) 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(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) 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 ! ! 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(KB_sta,KC_sta,KD_sta)
deallocate(Om1,X1,Y1,Om2,X2,Y2) 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 ! ! 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(KB_sta,KC_sta,KD_sta)
deallocate(Om1,X1,Y1,Om2,X2,Y2) deallocate(Om1,X1,Y1,Om2,X2,Y2)

View File

@ -76,8 +76,8 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
end do end do
KD_dyn(ij,kl) = 2d0*KD_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) = 2d0*ZD_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
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 do m=1,nS
num = (rho(i,k,m)*rho(j,l,m) - rho(j,k,m)*rho(i,l,m))/2 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) dem = - OmBSE - Om(m) + eGW(j) + eGW(l)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**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 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 end do
KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl) KD_dyn(ij,kl) = KD_dyn(ij,kl)
ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl) ZD_dyn(ij,kl) = ZD_dyn(ij,kl)
end do end do
end do end do