mirror of
https://github.com/pfloos/quack
synced 2025-01-03 01:55:57 +01:00
working on ufBSE
This commit is contained in:
Pierre-Francois Loos
parent
180ed63d74
commit
767db708c9
@ -256,10 +256,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy =',sum(EcBSE),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
|
||||
ERI,dipole_int,eT,eGT,EcBSE)
|
||||
subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
|
||||
ERI,dipole_int,eT,eGT,EcBSE)
|
||||
|
||||
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
|
||||
|
||||
@ -18,7 +18,7 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
logical,intent(in) :: triplet
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
@ -30,10 +30,10 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
integer,intent(in) :: nVVab
|
||||
integer,intent(in) :: nVVaa
|
||||
|
||||
double precision,intent(in) :: eT(nBas)
|
||||
double precision,intent(in) :: eGT(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
double precision,intent(in) :: eT(nOrb)
|
||||
double precision,intent(in) :: eGT(nOrb)
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||
|
||||
double precision,intent(in) :: Om1ab(nVVab)
|
||||
double precision,intent(in) :: X1ab(nVVab,nVVab)
|
||||
@ -41,16 +41,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
double precision,intent(in) :: Om2ab(nOOab)
|
||||
double precision,intent(in) :: X2ab(nVVab,nOOab)
|
||||
double precision,intent(in) :: Y2ab(nOOab,nOOab)
|
||||
double precision,intent(in) :: rho1ab(nBas,nBas,nVVab)
|
||||
double precision,intent(in) :: rho2ab(nBas,nBas,nOOab)
|
||||
double precision,intent(in) :: rho1ab(nOrb,nOrb,nVVab)
|
||||
double precision,intent(in) :: rho2ab(nOrb,nOrb,nOOab)
|
||||
double precision,intent(in) :: Om1aa(nVVaa)
|
||||
double precision,intent(in) :: X1aa(nVVaa,nVVaa)
|
||||
double precision,intent(in) :: Y1aa(nOOaa,nVVaa)
|
||||
double precision,intent(in) :: Om2aa(nOOaa)
|
||||
double precision,intent(in) :: X2aa(nVVaa,nOOaa)
|
||||
double precision,intent(in) :: Y2aa(nOOaa,nOOaa)
|
||||
double precision,intent(in) :: rho1aa(nBas,nBas,nVVaa)
|
||||
double precision,intent(in) :: rho2aa(nBas,nBas,nOOaa)
|
||||
double precision,intent(in) :: rho1aa(nOrb,nOrb,nVVaa)
|
||||
double precision,intent(in) :: rho2aa(nOrb,nOrb,nOOaa)
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -100,16 +100,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
|
||||
allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab))
|
||||
|
||||
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp)
|
||||
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp)
|
||||
|
||||
call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
call RGTpp_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab)
|
||||
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
|
||||
call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab)
|
||||
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
|
||||
|
||||
!----------------------------------------!
|
||||
! Compute T-matrix for alpha-alpha block !
|
||||
@ -120,16 +120,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
|
||||
allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa))
|
||||
|
||||
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp)
|
||||
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp)
|
||||
|
||||
call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
call RGTpp_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa)
|
||||
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
|
||||
call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa)
|
||||
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
|
||||
|
||||
!------------------!
|
||||
! Singlet manifold !
|
||||
@ -141,8 +141,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
|
||||
! Compute BSE singlet excitation energies
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
Aph(:,:) = Aph(:,:) + TAaa(:,:) + TAab(:,:)
|
||||
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) + TBab(:,:)
|
||||
@ -150,12 +150,12 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
|
||||
|
||||
call print_excitation_energies('phBSE@GTpp','singlet',nS,OmBSE)
|
||||
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
|
||||
call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
|
||||
|
||||
! Compute dynamic correction for BSE via renormalized perturbation theory
|
||||
|
||||
if(dBSE) &
|
||||
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, &
|
||||
dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
|
||||
|
||||
@ -171,8 +171,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
|
||||
! Compute BSE triplet excitation energies
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
Aph(:,:) = Aph(:,:) + TAaa(:,:) - TAab(:,:)
|
||||
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) - TBab(:,:)
|
||||
@ -180,12 +180,12 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
|
||||
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
|
||||
|
||||
call print_excitation_energies('phBSE@GTpp','triplet',nS,OmBSE)
|
||||
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
|
||||
call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
|
||||
|
||||
! Compute dynamic correction for BSE via renormalized perturbation theory
|
||||
|
||||
if(dBSE) &
|
||||
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
|
||||
Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, &
|
||||
dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
|
||||
end if
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
|
||||
! Unfold BSE@GW equations
|
||||
! Upfolded phBSE@GW (TDA singlets only!)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
@ -43,9 +43,9 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**********************************************'
|
||||
write(*,*)'| Unfolded BSE@GW calculation |'
|
||||
write(*,*)'**********************************************'
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)'* Upfolded phBSE@GW Calculation *'
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)
|
||||
|
||||
! TDA for W
|
||||
@ -71,11 +71,11 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
! Compute BSE supermatrix !
|
||||
!---------------------------!
|
||||
! !
|
||||
! | A -Ve -Vh | !
|
||||
! | A Jph Kph | !
|
||||
! | | !
|
||||
! H = | +Vh C2h2p 0 | !
|
||||
! H = | Kph C2h2p 0 | !
|
||||
! | | !
|
||||
! | +Ve 0 C2p2h | !
|
||||
! | Jph 0 C2p2h | !
|
||||
! !
|
||||
!---------------------------!
|
||||
|
||||
@ -102,21 +102,19 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Blocks Vp & Ve !
|
||||
!----------------!
|
||||
!------------------!
|
||||
! Blocks Jph & Kph !
|
||||
!------------------!
|
||||
|
||||
iajb=0
|
||||
ia = 0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nOrb-nR
|
||||
do j=nC+1,nO
|
||||
do b=nO+1,nOrb-nR
|
||||
iajb = iajb + 1
|
||||
ia = ia + 1
|
||||
|
||||
kc = 0
|
||||
do k=nC+1,nO
|
||||
do c=nO+1,nOrb-nR
|
||||
kc = kc + 1
|
||||
do k=nC+1,nO
|
||||
do c=nO+1,nOrn-nR
|
||||
do m=1,nS
|
||||
kcm = kcm + 1
|
||||
|
||||
tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i)
|
||||
H(n1h1p+iajb,kc ) = +tmp
|
||||
@ -56,25 +56,25 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
|
||||
dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
|
||||
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(k) - Om(m) + eGW(j)
|
||||
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(j)
|
||||
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(i)
|
||||
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)
|
||||
@ -106,25 +106,25 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
|
||||
kl = kl + 1
|
||||
|
||||
do m=1,nS
|
||||
dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
|
||||
dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
|
||||
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(k) - Om(m) + eGW(j)
|
||||
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(j)
|
||||
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(i)
|
||||
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)
|
||||
|
||||
227
src/GW/RGW_ppBSE_upfolded.f90
Normal file
227
src/GW/RGW_ppBSE_upfolded.f90
Normal file
@ -0,0 +1,227 @@
|
||||
subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
|
||||
|
||||
! Upfolded ppBSE@GW (TDA triplets only!)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: eHF(nOrb)
|
||||
double precision,intent(in) :: eGW(nOrb)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: s
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
integer :: ia,jb,kc,iajb,kcld
|
||||
integer,parameter :: maxH = 20
|
||||
double precision :: tmp
|
||||
|
||||
integer :: n1h1p,n2h2p,nH
|
||||
double precision,external :: Kronecker_delta
|
||||
|
||||
integer,allocatable :: order(:)
|
||||
double precision,allocatable :: H(:,:)
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**********************************************'
|
||||
write(*,*)'| Unfolded ppBSE@GW calculation |'
|
||||
write(*,*)'**********************************************'
|
||||
write(*,*)
|
||||
|
||||
! TDA for W
|
||||
|
||||
write(*,*) 'Tamm-Dancoff approximation by default!'
|
||||
write(*,*)
|
||||
|
||||
! Dimension of the supermatrix
|
||||
|
||||
n1h1p = nO*nV
|
||||
n2h2p = nO*nO*nV*nV
|
||||
nH = n1h1p + n2h2p + n2h2p
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH))
|
||||
|
||||
! Initialization
|
||||
|
||||
H(:,:) = 0d0
|
||||
|
||||
!----------------------------------------!
|
||||
! Compute BSE supermatrix !
|
||||
!----------------------------------------!
|
||||
! !
|
||||
! | D -M1 -M1 -M2 -M2 | !
|
||||
! | | !
|
||||
! | +M1 E1 0 0 0 | !
|
||||
! | | !
|
||||
! H = | +M1 0 E2 0 0 | !
|
||||
! | | !
|
||||
! | +M2 0 0 E3 0 | !
|
||||
! | | !
|
||||
! | +M2 0 0 0 E4 | !
|
||||
! !
|
||||
!----------------------------------------!
|
||||
|
||||
!---------!
|
||||
! Block D !
|
||||
!---------!
|
||||
|
||||
ij = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
|
||||
+ (ERI(i,j,k,l) - ERI(i,j,l,k))
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Blocks M1 !
|
||||
!----------------!
|
||||
|
||||
ijm = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
do m=1,nS
|
||||
ijm = ijm + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m)
|
||||
H(n2h+ijm,kl ) = +tmp
|
||||
H(kl ,n2h+ijm) = -tmp
|
||||
|
||||
tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m)
|
||||
H(n2h+1*n3h1p+ijm,kl ) = +tmp
|
||||
H(kl ,n2h+n3h1p+ijm) = -tmp
|
||||
|
||||
tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m)
|
||||
H(n2h+2*n3h1p+iajb,kc ) = +tmp
|
||||
H(kc ,n1h1p+iajb) = -tmp
|
||||
|
||||
tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m)
|
||||
H(n2h+3*n2h2p+iajb,kc ) = +tmp
|
||||
H(kc ,n1h1p+iajb) = -tmp
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------!
|
||||
! Block 2h2p !
|
||||
!-------------!
|
||||
|
||||
iajb = 0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nOrb-nR
|
||||
do j=nC+1,nO
|
||||
do b=nO+1,nOrb-nR
|
||||
iajb = iajb + 1
|
||||
|
||||
kcld = 0
|
||||
do k=nC+1,nO
|
||||
do c=nO+1,nOrb-nR
|
||||
do l=nC+1,nO
|
||||
do d=nO+1,nOrb-nR
|
||||
kcld = kcld + 1
|
||||
|
||||
tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
|
||||
+ 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d)
|
||||
H(n1h1p +iajb,n1h1p +kcld) = tmp
|
||||
H(n1h1p+n2h2p+iajb,n1h1p+n2h2p+kcld) = tmp
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------------------!
|
||||
! Diagonalize supermatrix !
|
||||
!-------------------------!
|
||||
|
||||
! call matout(nH,nH,H)
|
||||
|
||||
call diagonalize_general_matrix(nH,H,Om,X)
|
||||
|
||||
do s=1,nH
|
||||
order(s) = s
|
||||
end do
|
||||
|
||||
call quick_sort(Om,order,nH)
|
||||
call set_order(X,order,nH,nH)
|
||||
|
||||
!-----------------!
|
||||
! Compute weights !
|
||||
!-----------------!
|
||||
|
||||
Z(:) = 0d0
|
||||
do s=1,nH
|
||||
do ia=1,n1h1p
|
||||
Z(s) = Z(s) + X(ia,s)**2
|
||||
end do
|
||||
end do
|
||||
|
||||
!--------------!
|
||||
! Dump results !
|
||||
!--------------!
|
||||
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)' unfolded BSE excitation energies (eV) '
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
|
||||
'|','#','|','Omega (eV)','|','Z','|'
|
||||
write(*,*)'-------------------------------------------'
|
||||
|
||||
do s=1,min(nH,maxH)
|
||||
if(Z(s) > 1d-7) &
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',s,'|',Om(s)*HaToeV,'|',Z(s),'|'
|
||||
end do
|
||||
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
Loading…
Reference in New Issue
Block a user