LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

clean GTpp in progress but unfinished

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-10-01 16:14:33 +02:00
parent 75479f07f2
commit da6dada877
5 changed files with 158 additions and 252 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/GT/RG0T0pp.f90
View File

@ -273,7 +273,6 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
if(doppBSE) then if(doppBSE) then
call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, & call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI,dipole_int,eHF,eGT,EcBSE) ERI,dipole_int,eHF,eGT,EcBSE)
write(*,*) write(*,*)

116
src/GT/RGTpp_phBSE.f90
View File

@ -1,6 +1,5 @@
subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t,ERI,dipole_int,eT,eGT,EcBSE)
ERI,dipole_int,eT,eGT,EcBSE)
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
@ -25,39 +24,38 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS integer,intent(in) :: nS
integer,intent(in) :: nOOab integer,intent(in) :: nOOs
integer,intent(in) :: nOOaa integer,intent(in) :: nOOt
integer,intent(in) :: nVVab integer,intent(in) :: nVVs
integer,intent(in) :: nVVaa integer,intent(in) :: nVVt
double precision,intent(in) :: eT(nOrb) double precision,intent(in) :: eT(nOrb)
double precision,intent(in) :: eGT(nOrb) double precision,intent(in) :: eGT(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: Om1ab(nVVab) double precision,intent(in) :: Om1s(nVVs)
double precision,intent(in) :: X1ab(nVVab,nVVab) double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1ab(nOOab,nVVab) double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: Om2ab(nOOab) double precision,intent(in) :: Om2s(nOOs)
double precision,intent(in) :: X2ab(nVVab,nOOab) double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2ab(nOOab,nOOab) double precision,intent(in) :: Y2s(nOOs,nOOs)
double precision,intent(in) :: rho1ab(nOrb,nOrb,nVVab) double precision,intent(in) :: rho1s(nOrb,nOrb,nVVs)
double precision,intent(in) :: rho2ab(nOrb,nOrb,nOOab) double precision,intent(in) :: rho2s(nOrb,nOrb,nOOs)
double precision,intent(in) :: Om1aa(nVVaa) double precision,intent(in) :: Om1t(nVVt)
double precision,intent(in) :: X1aa(nVVaa,nVVaa) double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1aa(nOOaa,nVVaa) double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: Om2aa(nOOaa) double precision,intent(in) :: Om2t(nOOt)
double precision,intent(in) :: X2aa(nVVaa,nOOaa) double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2aa(nOOaa,nOOaa) double precision,intent(in) :: Y2t(nOOt,nOOt)
double precision,intent(in) :: rho1aa(nOrb,nOrb,nVVaa) double precision,intent(in) :: rho1t(nOrb,nOrb,nVVt)
double precision,intent(in) :: rho2aa(nOrb,nOrb,nOOaa) double precision,intent(in) :: rho2t(nOrb,nOrb,nOOt)
! Local variables ! Local variables
logical :: dRPA logical :: dRPA = .false.
integer :: ispin integer :: ispin
integer :: iblock
double precision,allocatable :: Bpp(:,:) double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Cpp(:,:) double precision,allocatable :: Cpp(:,:)
@ -67,8 +65,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
double precision,allocatable :: Bph(:,:) double precision,allocatable :: Bph(:,:)
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision,allocatable :: TAab(:,:),TBab(:,:) double precision,allocatable :: TAs(:,:),TBs(:,:)
double precision,allocatable :: TAaa(:,:),TBaa(:,:) double precision,allocatable :: TAt(:,:),TBt(:,:)
double precision,allocatable :: OmBSE(:) double precision,allocatable :: OmBSE(:)
double precision,allocatable :: XpY_BSE(:,:) double precision,allocatable :: XpY_BSE(:,:)
double precision,allocatable :: XmY_BSE(:,:) double precision,allocatable :: XmY_BSE(:,:)
@ -79,7 +77,7 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
! Memory allocation ! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),TAab(nS,nS),TBab(nS,nS),TAaa(nS,nS),TBaa(nS,nS), & allocate(Aph(nS,nS),Bph(nS,nS),TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), &
OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
!-----! !-----!
@ -91,45 +89,43 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
write(*,*) write(*,*)
end if end if
!---------------------------------------! !------------------------------------!
! Compute T-matrix for alpha-beta block ! ! Compute T-matrix for singlet block !
!---------------------------------------! !------------------------------------!
ispin = 1 ispin = 1
iblock = 3
allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab) call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Om1s,rho1s,Om2s,rho2s,TAs)
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab) if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Om1s,rho1s,Om2s,rho2s,TBs)
!----------------------------------------! !------------------------------------!
! Compute T-matrix for alpha-alpha block ! ! Compute T-matrix for triplet block !
!----------------------------------------! !------------------------------------!
ispin = 2 ispin = 2
iblock = 4
allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa) call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Om1t,rho1t,Om2t,rho2t,TAt)
if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa) if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Om1t,rho1t,Om2t,rho2t,TBt)
!------------------! !------------------!
! Singlet manifold ! ! Singlet manifold !
@ -144,8 +140,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) 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) if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
Aph(:,:) = Aph(:,:) + TAaa(:,:) + TAab(:,:) Aph(:,:) = Aph(:,:) + TAt(:,:) ! TAt(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) + TBab(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) ! TBt(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
@ -155,9 +151,9 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
! Compute dynamic correction for BSE via renormalized perturbation theory ! Compute dynamic correction for BSE via renormalized perturbation theory
if(dBSE) & if(dBSE) &
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & Om1s,Om2s,Om1t,Om2t,rho1s,rho2s,rho1t,rho2t,eT,eGT, &
dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAs,TAt)
end if end if
@ -174,8 +170,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) 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) if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
Aph(:,:) = Aph(:,:) + TAaa(:,:) - TAab(:,:) Aph(:,:) = Aph(:,:) + 1d0*TAt(:,:) - TAs(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) - TBab(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + 1d0*TBt(:,:) - TBs(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
@ -185,9 +181,9 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
! Compute dynamic correction for BSE via renormalized perturbation theory ! Compute dynamic correction for BSE via renormalized perturbation theory
if(dBSE) & if(dBSE) &
call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & Om1s,Om2s,Om1t,Om2t,rho1s,rho2s,rho1t,rho2t,eT,eGT, &
dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAs,TAt)
end if end if
if(exchange_kernel) then if(exchange_kernel) then

223
src/GT/RGTpp_ppBSE.f90
View File

@ -1,5 +1,4 @@
subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, &
Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
ERI,dipole_int,eT,eGT,EcBSE) ERI,dipole_int,eT,eGT,EcBSE)
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
@ -23,59 +22,90 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOOab integer,intent(in) :: nOOs
integer,intent(in) :: nOOaa integer,intent(in) :: nOOt
integer,intent(in) :: nVVab integer,intent(in) :: nVVs
integer,intent(in) :: nVVaa integer,intent(in) :: nVVt
double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas) double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Om1ab(nVVab)
double precision,intent(in) :: X1ab(nVVab,nVVab)
double precision,intent(in) :: Y1ab(nOOab,nVVab)
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) :: 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)
! Local variables ! Local variables
integer :: ispin integer :: ispin
integer :: iblock
integer :: nOOs
integer :: nOOt
integer :: nVVs
integer :: nVVt
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:) double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
double precision,allocatable :: TBab(:,:),TCab(:,:),TDab(:,:)
double precision,allocatable :: TBaa(:,:),TCaa(:,:),TDaa(:,:)
double precision,allocatable :: Om1s(:),Om1t(:) double precision,allocatable :: Om1s(:),Om1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Om2s(:),Om2t(:) double precision,allocatable :: Om2s(:),Om2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:) double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
double precision,allocatable :: TBs(:,:),TCs(:,:),TDs(:,:)
double precision,allocatable :: TBt(:,:),TCt(:,:),TDt(:,:)
! Output variables ! Output variables
double precision,intent(out) :: EcBSE(nspin) double precision,intent(out) :: EcBSE(nspin)
!------------------------------------!
! Compute T-matrix for singlet block !
!------------------------------------!
ispin = 1
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBs(nVVs,nOOs),TCs(nVVs,nVVs),TDs(nOOs,nOOs))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, &
Om1s,rho1s,Om2s,rho2s,TBs)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, &
Om1s,rho1s,Om2s,rho2s,TCs)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, &
Om1s,rho1s,Om2s,rho2s,TDs)
deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s)
!------------------------------------!
! Compute T-matrix for triplet block !
!------------------------------------!
ispin = 2
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBt(nVVt,nOOt),TCt(nVVt,nVVt),TDt(nOOt,nOOt))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, &
Om1t,rho1t,Om2t,rho2t,TBt)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, &
Om1t,rho1t,Om2t,rho2t,TCt)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, &
Om1t,rho1t,Om2t,rho2t,TDt)
deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t)
!------------------! !------------------!
! Singlet manifold ! ! Singlet manifold !
@ -85,63 +115,6 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
ispin = 1 ispin = 1
nOOs = nO*(nO+1)/2
nVVs = nV*(nV+1)/2
!---------------------------------------!
! Compute T-matrix for alpha-beta block !
!---------------------------------------!
iblock = 3
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)
call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBab(nVVs,nOOs),TCab(nVVs,nVVs),TDab(nOOs,nOOs))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TBab)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TCab)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TDab)
!----------------------------------------!
! Compute T-matrix for alpha-alpha block !
!----------------------------------------!
iblock = 4
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)
call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBaa(nVVs,nOOs),TCaa(nVVs,nVVs),TDaa(nOOs,nOOs))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
!----------------------------------!
! pp/hh sectors for singlet states !
!----------------------------------!
EcBSE(ispin) = 0d0
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
@ -149,15 +122,15 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) - TBab(:,:) - TBaa(:,:) Bpp(:,:) = Bpp(:,:) - TBs(:,:) - TBt(:,:)
Cpp(:,:) = Cpp(:,:) - TCab(:,:) - TCaa(:,:) Cpp(:,:) = Cpp(:,:) - TCs(:,:) - TCt(:,:)
Dpp(:,:) = Dpp(:,:) - TDab(:,:) - TDaa(:,:) Dpp(:,:) = Dpp(:,:) - TDs(:,:) - TDt(:,:)
call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin)) call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s) call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,TBab,TCab,TDab,TBaa,TCaa,TDaa,Bpp,Cpp,Dpp) deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,Bpp,Cpp,Dpp)
end if end if
@ -169,80 +142,24 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
ispin = 2 ispin = 2
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
!---------------------------------------!
! Compute T-matrix for alpha-beta block !
!---------------------------------------!
iblock = 3
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)
call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBab(nVVt,nOOt),TCab(nVVt,nVVt),TDab(nOOt,nOOt))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TBab)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TCab)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0, &
Om1ab,rho1ab,Om2ab,rho2ab,TDab)
!----------------------------------------!
! Compute T-matrix for alpha-alpha block !
!----------------------------------------!
iblock = 4
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)
call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp)
allocate(TBaa(nVVt,nOOt),TCaa(nVVt,nVVt),TDaa(nOOt,nOOt))
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0, &
Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
!----------------------------------!
! pp/hh sectors for triplet states !
!----------------------------------!
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
Bpp(:,:) = Bpp(:,:) + TBab(:,:) - TBaa(:,:) Bpp(:,:) = Bpp(:,:) + TBs(:,:) - TBt(:,:)
Cpp(:,:) = Cpp(:,:) + TCab(:,:) - TCaa(:,:) Cpp(:,:) = Cpp(:,:) + TCs(:,:) - TCt(:,:)
Dpp(:,:) = Dpp(:,:) + TDab(:,:) - TDaa(:,:) Dpp(:,:) = Dpp(:,:) + TDs(:,:) - TDt(:,:)
call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin)) call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t) call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t,TBab,TCab,TDab,TBaa,TCaa,TDaa,Bpp,Cpp,Dpp) deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t,Bpp,Cpp,Dpp)
end if end if

29
src/GT/evRGTpp.f90
View File

@ -48,7 +48,6 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
double precision :: rcond double precision :: rcond
double precision :: Conv double precision :: Conv
integer :: ispin integer :: ispin
integer :: iblock
integer :: nOOs,nOOt integer :: nOOs,nOOt
integer :: nVVs,nVVt integer :: nVVs,nVVt
double precision :: EcGM double precision :: EcGM
@ -91,8 +90,8 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
! Dimensions of the pp-RPA linear reponse matrices ! Dimensions of the pp-RPA linear reponse matrices
nOOs = nO*nO nOOs = nO*(nO+1)/2
nVVs = nV*nV nVVs = nV*(nV+1)/2
nOOt = nO*(nO-1)/2 nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2 nVVt = nV*(nV-1)/2
@ -131,15 +130,14 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
!---------------------------------------------- !----------------------------------------------
ispin = 1 ispin = 1
iblock = 3
! Compute linear response ! Compute linear response
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
@ -150,32 +148,31 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
!---------------------------------------------- !----------------------------------------------
ispin = 2 ispin = 2
iblock = 4
! Compute linear response ! Compute linear response
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = 1d0*EcRPA(1)
EcRPA(2) = 3d0*EcRPA(2) EcRPA(2) = 3d0*EcRPA(2)
!---------------------------------------------- !----------------------------------------------
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy
!---------------------------------------------- !----------------------------------------------
iblock = 3 ispin = 3
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) call RGTpp_excitation_density(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 4 ispin = 4
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) call RGTpp_excitation_density(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
!---------------------------------------------- !----------------------------------------------
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy

29
src/GT/qsRGTpp.f90
View File

@ -53,7 +53,6 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
integer :: nSCF integer :: nSCF
integer :: nBas_Sq integer :: nBas_Sq
integer :: ispin integer :: ispin
integer :: iblock
integer :: n_diis integer :: n_diis
double precision :: ET double precision :: ET
double precision :: EV double precision :: EV
@ -108,8 +107,8 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Dimensions of the pp-RPA linear reponse matrices ! Dimensions of the pp-RPA linear reponse matrices
nOOs = nO*nO nOOs = nO*(nO+1)/2
nVVs = nV*nV nVVs = nV*(nV+1)/2
nOOt = nO*(nO-1)/2 nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2 nVVt = nV*(nV-1)/2
@ -196,41 +195,39 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute linear response ! Compute linear response
ispin = 1 ispin = 1
iblock = 3
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
ispin = 2 ispin = 2
iblock = 4
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) if(.not.TDA_T) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = 1d0*EcRPA(1)
EcRPA(2) = 3d0*EcRPA(2) EcRPA(2) = 3d0*EcRPA(2)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
iblock = 3 ispin = 3
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) call RGTpp_excitation_density(ispin,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 4 ispin = 4
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) call RGTpp_excitation_density(ispin,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
if(regularize) then if(regularize) then
call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)