LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 10:05:59 +01:00
Code Releases Activity

ppBSE GT

Browse Source
This commit is contained in:
Antoine Marie 2024-10-30 09:36:36 +01:00
parent 61288e474d
commit 6c37311e38
25 changed files with 771 additions and 462 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
src/GT/GG0T0pp.f90
View File

@ -1,4 +1,4 @@
subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE, & subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -13,7 +13,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
logical,intent(in) :: dophBSE logical,intent(in) :: dophBSE,dophBSE2
logical,intent(in) :: doppBSE logical,intent(in) :: doppBSE
logical,intent(in) :: TDA_T logical,intent(in) :: TDA_T
logical,intent(in) :: TDA logical,intent(in) :: TDA
@ -37,7 +37,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Local variables ! Local variables
logical :: print_T = .false. logical :: print_T = .true.
integer :: nOO integer :: nOO
integer :: nVV integer :: nVV
@ -72,7 +72,6 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
write(*,*)'**********************************' write(*,*)'**********************************'
write(*,*) write(*,*)
! TDA for T ! TDA for T
if(TDA_T) then if(TDA_T) then
@ -104,7 +103,9 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Compute linear response ! Compute linear response
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
Bpp(:,:) = 0d0
Cpp(:,:) = 0d0
Dpp(:,:) = 0d0
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
@ -114,12 +115,11 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
if(print_T) call print_excitation_energies('ppRPA@GHF','2p',nVV,Om1) if(print_T) call print_excitation_energies('ppRPA@GHF','2p',nVV,Om1)
if(print_T) call print_excitation_energies('ppRPA@FHF','2h',nOO,Om2) if(print_T) call print_excitation_energies('ppRPA@GHF','2h',nOO,Om2)
!---------------------------------------------- !----------------------------------------------
! Compute excitation densities ! Compute excitation densities
!---------------------------------------------- !----------------------------------------------
call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
!---------------------------------------------- !----------------------------------------------
@ -170,6 +170,24 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA) call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA)
!----------------------------------------------
! ppBSE calculation
!----------------------------------------------
if(doppBSE) then
call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
ERI,dipole_int,eHF,eGT,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy = ',EcBSE,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF total energy = ',ENuc + EGHF + EcBSE,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
! Testing zone ! Testing zone
if(dotest) then if(dotest) then

117
src/GT/GGT.f90 Normal file
View File

@ -0,0 +1,117 @@
subroutine GGT(dotest,doG0T0,doevGT,doqsGT,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
TDA_T,TDA,dBSE,dTDA,linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! GT module
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doG0T0
logical,intent(in) :: doevGT
logical,intent(in) :: doqsGT
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: TDA_T
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: doppBSE
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: doSRG
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas2
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas2)
double precision,intent(in) :: cHF(nBas2,nBas2)
double precision,intent(in) :: PHF(nBas2,nBas2)
double precision,intent(in) :: S(nBas2,nBas2)
double precision,intent(in) :: T(nBas2,nBas2)
double precision,intent(in) :: V(nBas2,nBas2)
double precision,intent(in) :: Hc(nBas2,nBas2)
double precision,intent(in) :: X(nBas2,nBas2)
double precision,intent(in) :: ERI_AO(nBas2,nBas2,nBas2,nBas2)
double precision,intent(in) :: ERI(nBas2,nBas2,nBas2,nBas2)
double precision,intent(in) :: dipole_int_AO(nBas2,nBas2,ncart)
double precision,intent(in) :: dipole_int(nBas2,nBas2,ncart)
! Local variables
double precision :: start_GT ,end_GT ,t_GT
!------------------------------------------------------------------------
! Perform G0T0 calculatiom
!------------------------------------------------------------------------
if(doG0T0) then
call wall_time(start_GT)
call GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGT calculation
!------------------------------------------------------------------------
! if(doevGT) then
! call wall_time(start_GT)
! call evGGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
! linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! call wall_time(end_GT)
! t_GT = end_GT - start_GT
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds'
! write(*,*)
! end if
!------------------------------------------------------------------------
! Perform qsGT calculation
!------------------------------------------------------------------------
! if(doqsGT) then
! call wall_time(start_GT)
! call qsGGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
! eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI, &
! dipole_int_AO,dipole_int,PHF,cHF,eHF)
! call wall_time(end_GT)
! t_GT = end_GT - start_GT
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_GT,' seconds'
! write(*,*)
! end if
end subroutine

74
src/GT/GGT_Tmatrix.f90 Normal file
View File

@ -0,0 +1,74 @@
subroutine GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T)
! Compute the T-matrix tensor elements
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eGF(nOrb)
double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nOrb,nOrb,nOO)
! Local variables
integer :: p,q,r,s
integer :: c,d,k,l
integer :: kl,cd
! Output variables
double precision,intent(out) :: T(nOrb,nOrb,nOrb,nOrb)
! Initialization
T(:,:,:,:) = 0d0
! Start by building the tensor elements of T
! This is probabbly not a good idea because this tensor is really large
!$OMP PARALLEL &
!$OMP SHARED(nC,nO,nOrb,nR,T,ERI,rho1,rho2,Om1,Om2) &
!$OMP PRIVATE(p,q,r,s,c,d,cd,k,l,kl) &
!$OMP DEFAULT(NONE)
!$OMP DO
do s=nC+1,nOrb-nR
do r=nC+1,nOrb-nR
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r)
cd = 0
do c = nO+1, nOrb-nR
do d = c+1, nOrb-nR
cd = cd + 1
T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd) * rho1(r,s,cd) / Om1(cd)
end do ! d
end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
kl = kl + 1
T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl) * rho2(r,s,kl) / Om2(kl)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end subroutine

3
src/GT/GGTpp_excitation_density.f90
View File

@ -1,4 +1,4 @@
subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) subroutine GGTpp_excitation_density(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities for T-matrix self-energy ! Compute excitation densities for T-matrix self-energy
@ -12,7 +12,6 @@ subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
! Input variables ! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO

127
src/GT/GGTpp_ppBSE.f90 Normal file
View File

@ -0,0 +1,127 @@
subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
ERI,dipole_int,eT,eGT,EcBSE)
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: TDA_T
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
double precision,intent(in) :: eta
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
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)
! Local variables
double precision :: EcRPA(nspin)
double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
double precision,allocatable :: Om1(:), Om2(:)
double precision,allocatable :: X1(:,:), X2(:,:)
double precision,allocatable :: Y1(:,:), Y2(:,:)
double precision,allocatable :: rho1(:,:,:), rho2(:,:,:)
double precision,allocatable :: KB_sta(:,:),KC_sta(:,:),KD_sta(:,:)
double precision,allocatable :: T(:,:,:,:)
! Output variables
double precision,intent(out) :: EcBSE
!----------------------------------------------
! Compute linear response
!----------------------------------------------
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
Bpp(:,:) = 0d0
Cpp(:,:) = 0d0
Dpp(:,:) = 0d0
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eT,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
deallocate(Bpp,Cpp,Dpp)
!----------------------------------------------
! Compute excitation densities
!----------------------------------------------
allocate(rho1(nOrb,nOrb,nVV),rho2(nOrb,nOrb,nOO))
call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
deallocate(X1,Y1,X2,Y2)
!----------------------------------------------
! Compute new ppRPA block
!----------------------------------------------
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
!----------------------------------------------
! Compute T matrix tensor
!----------------------------------------------
allocate(T(nOrb,nOrb,nOrb,nOrb))
call GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T)
!----------------------------------------------
! Compute kernels
!----------------------------------------------
allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KC_sta)
call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KD_sta)
if(.not.TDA_T) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KB_sta)
deallocate(Om1,Om2,rho1,rho2)
! Deallocate the 4-tensor T
deallocate(T)
!----------------------------------------------
! Diagonalize ppBSE
!----------------------------------------------
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
!----------------------------------------------------!
! Compute the dynamical screening at the ppBSE level !
!----------------------------------------------------!
! TODO
deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
end subroutine

56
src/GT/GGTpp_ppBSE_static_kernel_B.f90
View File

@ -1,42 +1,48 @@
subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB) subroutine GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KB_sta)
! Compute the VVOO block of the static T-matrix ! Compute the VVOO of the T-matrix static pp kernel
implicit none implicit none
include 'parameters.h' include 'parameters.h'
! Input variables ! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nOrb
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eGF(nOrb)
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: rho2(nOrb,nOrb,nOO)
double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb)
! Local variables ! Local variables
double precision :: chi double precision :: dem,num
double precision :: eps integer :: p,q,r,s,e,m
integer :: i,j,a,b,ij,ab,cd,kl integer :: i,j,a,b
integer :: ij,ab,kl,cd
! Output variables ! Output variables
double precision,intent(out) :: TB(nVVx,nOOx) double precision,intent(out) :: KB_sta(nVV,nOO)
! Initialization
KB_sta(:,:) = 0d0
! Computing the kernel
! This is the same code as for the GF(2) kernel with elements T instead of ERI
ab = 0 ab = 0
do a=nO+1,nBas-nR do a=nO+1,nOrb-nR
do b=a+1,nBas-nR do b=a+1,nOrb-nR
ab = ab + 1 ab = ab + 1
ij = 0 ij = 0
@ -44,19 +50,15 @@ subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
do j=i+1,nO do j=i+1,nO
ij = ij + 1 ij = ij + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nOrb-nR
do cd=1,nVV dem = eGF(m) - eGF(e)
eps = + Om1(cd) num = (T(a,m,i,e) - T(a,m,e,i)) * (T(e,b,m,j) - T(e,b,j,m))
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) num = num + (T(a,e,i,m) - T(a,e,m,i)) * (T(m,b,e,j) - T(m,b,j,e))
end do num = num - (T(b,m,i,e) - T(b,m,e,i)) * (T(e,a,m,j) - T(e,a,j,m))
num = num - (T(b,e,i,m) - T(b,e,m,i)) * (T(m,a,e,j) - T(m,a,j,e))
do kl=1,nOO KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
eps = - Om2(kl)
chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2)
end do
TB(ab,ij) = lambda*chi
end do end do
end do end do
@ -64,5 +66,7 @@ subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
end do end do
end do end do
end do
end do
end subroutine end subroutine

60
src/GT/GGTpp_ppBSE_static_kernel_C.f90
View File

@ -1,64 +1,70 @@
subroutine GGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC) subroutine GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KC_sta)
! Compute the VVVV block of the static T-matrix ! Compute the VVVV block of the T-matrix static pp kernels
implicit none implicit none
include 'parameters.h' include 'parameters.h'
! Input variables ! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nOrb
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eGF(nOrb)
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: rho2(nOrb,nOrb,nOO)
double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb)
! Local variables ! Local variables
double precision,external :: Kronecker_delta double precision :: dem,num
double precision :: chi integer :: p,q,r,s,e,m
double precision :: eps integer :: a,b,c,d,k,l
integer :: a,b,c,d,ab,cd,ef,mn integer :: ab,kl,cd
! Output variables ! Output variables
double precision,intent(out) :: TC(nVVx,nVVx) double precision,intent(out) :: KC_sta(nVV,nVV)
! Initialization
KC_sta(:,:) = 0d0
! Computing the kernel
! This is the same code as for the GF(2) kernel with elements T instead of ERI
ab = 0 ab = 0
do a=nO+1,nBas-nR do a=nO+1,nOrb-nR
do b=a+1,nBas-nR do b=a+1,nOrb-nR
ab = ab + 1 ab = ab + 1
cd = 0 cd = 0
do c=nO+1,nBas-nR do c=nO+1,nOrb-nR
do d=c+1,nBas-nR do d=c+1,nOrb-nR
cd = cd + 1 cd = cd + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nOrb-nR
dem = eGF(m) - eGF(e)
num = (T(a,m,c,e) - T(a,m,e,c)) * (T(e,b,m,d) - T(e,b,d,m))
num = num + (T(a,e,c,m) - T(a,e,m,c)) * (T(m,b,e,d) - T(m,b,d,e))
num = num - (T(b,m,c,e) - T(b,m,e,c)) * (T(e,a,m,d) - T(e,a,d,m))
num = num - (T(b,e,c,m) - T(b,e,m,c)) * (T(m,a,e,d) - T(m,a,d,e))
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
do ef=1,nVV
eps = + Om1(ef)
chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2)
end do end do
do mn=1,nOO
eps = - Om2(mn)
chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2)
end do end do
TC(ab,cd) = lambda*chi
end do end do
end do end do

50
src/GT/GGTpp_ppBSE_static_kernel_D.f90
View File

@ -1,39 +1,45 @@
subroutine GGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD) subroutine GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KD_sta)
! Compute the OOOO block of the static T-matrix ! Compute the OOOO block of the T-matrix static pp kernel
implicit none implicit none
include 'parameters.h' include 'parameters.h'
! Input variables ! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nOrb
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eGF(nOrb)
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: rho2(nOrb,nOrb,nOO)
double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb)
! Local variables ! Local variables
double precision :: chi double precision :: dem,num
double precision :: eps integer :: p,q,r,s,e,m
integer :: i,j,k,l,ij,kl,ef,mn integer :: i,j,k,l
integer :: ij,kl,cd
! Output variables ! Output variables
double precision,intent(out) :: TD(nOOx,nOOx) double precision,intent(out) :: KD_sta(nOO,nOO)
! Initialization
KD_sta(:,:) = 0d0
! Computing the kernel
! This is the same code as for the GF(2) kernel with elements T instead of ERI
ij = 0 ij = 0
do i=nC+1,nO do i=nC+1,nO
do j=i+1,nO do j=i+1,nO
@ -44,24 +50,20 @@ subroutine GGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
do l=k+1,nO do l=k+1,nO
kl = kl + 1 kl = kl + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nOrb-nR
do ef=1,nVV dem = eGF(m) - eGF(e)
eps = + Om1(ef) num = T(i,m,k,e) * T(e,j,m,l) + T(i,e,k,m) * T(m,j,e,l)
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2) num = num - T(j,m,k,e) * T(e,i,m,l) - T(j,e,k,m) * T(m,i,e,l)
end do
do mn=1,nOO
eps = - Om2(mn)
chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2)
end do
TD(ij,kl) = lambda*chi
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
end do end do
end do end do
end do end do
end do end do
end do
end do
end subroutine end subroutine

2
src/GT/RG0T0pp.f90
View File

@ -212,6 +212,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR_D(isp_T,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(isp_T)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
@ -224,6 +225,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) call ppLR_D(isp_T,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(isp_T)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)

3
src/GT/RGT.f90
View File

@ -1,6 +1,3 @@
! ---
subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, & doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, &

8
src/GT/RGTpp_SigC.f90
View File

@ -42,12 +42,12 @@ double precision function RGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV
do cd=1,nVVs do cd=1,nVVs
eps = w + e(i) - Om1s(cd) eps = w + e(i) - Om1s(cd)
RGTpp_SigC = RGTpp_SigC + rho1s(p,i,cd)**2*eps/(eps**2 + eta**2) RGTpp_SigC = RGTpp_SigC + (1d0/2d0)*rho1s(p,i,cd)**2*eps/(eps**2 + eta**2)
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = w + e(i) - Om1t(cd) eps = w + e(i) - Om1t(cd)
RGTpp_SigC = RGTpp_SigC + rho1t(p,i,cd)**2*eps/(eps**2 + eta**2) RGTpp_SigC = RGTpp_SigC + (3d0/2d0)*rho1t(p,i,cd)**2*eps/(eps**2 + eta**2)
end do end do
end do end do
@ -60,12 +60,12 @@ double precision function RGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV
do kl=1,nOOs do kl=1,nOOs
eps = w + e(a) - Om2s(kl) eps = w + e(a) - Om2s(kl)
RGTpp_SigC = RGTpp_SigC + rho2s(p,a,kl)**2*eps/(eps**2 + eta**2) RGTpp_SigC = RGTpp_SigC + (1d0/2d0)*rho2s(p,a,kl)**2*eps/(eps**2 + eta**2)
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = w + e(a) - Om2t(kl) eps = w + e(a) - Om2t(kl)
RGTpp_SigC = RGTpp_SigC + rho2t(p,a,kl)**2*eps/(eps**2 + eta**2) RGTpp_SigC = RGTpp_SigC + (3d0/2d0)*rho2t(p,a,kl)**2*eps/(eps**2 + eta**2)
end do end do
end do end do

9
src/GT/RGTpp_dSigC.f90
View File

@ -42,31 +42,30 @@ double precision function RGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do cd=1,nVVs do cd=1,nVVs
eps = w + e(i) - Om1s(cd) eps = w + e(i) - Om1s(cd)
RGTpp_dSigC = RGTpp_dSigC - rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 RGTpp_dSigC = RGTpp_dSigC - (1d0/2d0)*rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = w + e(i) - Om1t(cd) eps = w + e(i) - Om1t(cd)
RGTpp_dSigC = RGTpp_dSigC - rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 RGTpp_dSigC = RGTpp_dSigC - (3d0/2d0)*rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
!---------------------------------------------- !----------------------------------------------
! Virtual part of the T-matrix self-energy ! Virtual part of the T-matrix self-energy
!
!---------------------------------------------- !----------------------------------------------
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do kl=1,nOOs do kl=1,nOOs
eps = w + e(a) - Om2s(kl) eps = w + e(a) - Om2s(kl)
RGTpp_dSigC = RGTpp_dSigC - rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 RGTpp_dSigC = RGTpp_dSigC - (1d0/2d0)*rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = w + e(a) - Om2t(kl) eps = w + e(a) - Om2t(kl)
RGTpp_dSigC = RGTpp_dSigC - rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 RGTpp_dSigC = RGTpp_dSigC - (3d0/2d0)*rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do

20
src/GT/RGTpp_excitation_density.f90
View File

@ -66,7 +66,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
cd = cd + 1 cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) & rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/ & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/ &
sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) sqrt(1d0 + Kronecker_delta(c,d))
end do end do
end do end do
@ -76,7 +76,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
kl = kl + 1 kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) & rho1(p,q,ab) = rho1(p,q,ab) &
+ (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/ & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/ &
sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) sqrt(1d0 + Kronecker_delta(k,l))
end do end do
end do end do
@ -93,7 +93,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
cd = cd + 1 cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) & rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/ & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/ &
sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) sqrt(1d0 + Kronecker_delta(c,d))
end do end do
end do end do
@ -103,7 +103,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
kl = kl + 1 kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) & rho2(p,q,ij) = rho2(p,q,ij) &
+ (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/ & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/ &
sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) sqrt(1d0 + Kronecker_delta(k,l))
end do end do
end do end do
@ -148,7 +148,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
cd = cd + 1 cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) & rho1(p,q,ab) = rho1(p,q,ab) &
+ sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)
end do ! d end do ! d
end do ! c end do ! c
@ -159,7 +159,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
kl = kl + 1 kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) & rho1(p,q,ab) = rho1(p,q,ab) &
+ sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab)
end do ! l end do ! l
end do ! k end do ! k
end do ! b end do ! b
@ -179,7 +179,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
cd = cd + 1 cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) & rho2(p,q,ij) = rho2(p,q,ij) &
+ sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij)
end do ! d end do ! d
end do ! c end do ! c
@ -190,7 +190,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
kl = kl + 1 kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) & rho2(p,q,ij) = rho2(p,q,ij) &
+ sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)
end do ! l end do ! l
end do ! k end do ! k
end do ! j end do ! j
@ -249,8 +249,8 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
deallocate(ERI_1, ERI_2) deallocate(ERI_1, ERI_2)
rho1 = rho1*sqrt(3d0/2d0) rho1 = rho1
rho2 = rho2*sqrt(3d0/2d0) rho2 = rho2
endif endif
endif endif

144
src/GT/RGTpp_ppBSE.f90
View File

@ -34,7 +34,7 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
! Local variables ! Local variables
integer :: ispin integer :: ispin, isp_T
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:) double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
@ -46,66 +46,74 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
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 :: rho2s(:,:,:),rho2t(:,:,:)
double precision,allocatable :: TBs(:,:),TCs(:,:),TDs(:,:) double precision,allocatable :: KB_sta(:,:)
double precision,allocatable :: TBt(:,:),TCt(:,:),TDt(:,:) double precision,allocatable :: KC_sta(:,:)
double precision,allocatable :: KD_sta(:,:)
double precision,allocatable :: Taaaa(:,:,:,:),Tabab(:,:,:,:),Tbaab(:,:,:,:)
! Output variables ! Output variables
double precision,intent(out) :: EcBSE(nspin) double precision,intent(out) :: EcBSE(nspin)
!------------------------------------! !---------------------------------
! Compute T-matrix for singlet block ! ! Compute ppRPA excitation density
!------------------------------------! !---------------------------------
ispin = 1 isp_T = 1
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs))
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp) allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
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(isp_T))
! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
deallocate(Bpp,Cpp,Dpp) allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs))
allocate(TBs(nVVs,nOOs),TCs(nVVs,nVVs),TDs(nOOs,nOOs)) call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, & deallocate(X1s,Y1s,X2s,Y2s,Bpp,Cpp,Dpp)
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) isp_T = 2
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt))
allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
!------------------------------------! if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
! Compute T-matrix for triplet block ! call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
!------------------------------------! call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
ispin = 2 call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & allocate(rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) deallocate(X1t,Y1t,X2t,Y2t,Bpp,Cpp,Dpp)
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)) !---------------------------------
! Compute T matrix elements
!---------------------------------
deallocate(Bpp,Cpp,Dpp) ! This correspond to the alpha-alpha-alpha-alpha elements
allocate(TBt(nVVt,nOOt),TCt(nVVt,nVVt),TDt(nOOt,nOOt)) isp_T = 1
allocate(Taaaa(nBas,nBas,nBas,nBas))
call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Taaaa)
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, & ! This correspond to the alpha-alpha-alpha-alpha elements
Om1t,rho1t,Om2t,rho2t,TBt) isp_T = 2
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, & allocate(Tabab(nBas,nBas,nBas,nBas))
Om1t,rho1t,Om2t,rho2t,TCt) call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tabab)
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) ! This correspond to the alpha-alpha-alpha-alpha elements
isp_T = 3
allocate(Tbaab(nBas,nBas,nBas,nBas))
call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tbaab)
deallocate(Om1s,Om2s,Om1t,Om2t,rho1s,rho2s,rho1t,rho2t)
!------------------! !------------------!
! Singlet manifold ! ! Singlet manifold !
@ -115,22 +123,34 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
ispin = 1 ispin = 1
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))
Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
! Compute BSE excitation energies
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
allocate(KB_sta(nVVs,nOOs),KC_sta(nVVs,nVVs),KD_sta(nOOs,nOOs))
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(:,:) - TBs(:,:) - TBt(:,:) if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, &
Cpp(:,:) = Cpp(:,:) - TCs(:,:) - TCt(:,:) Taaaa,Tabab,Tbaab,KB_sta)
Dpp(:,:) = Dpp(:,:) - TDs(:,:) - TDt(:,:) call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, &
Taaaa,Tabab,Tbaab,KC_sta)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, &
Taaaa,Tabab,Tbaab,KD_sta)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
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,Bpp,Cpp,Dpp) deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
end if end if
@ -144,18 +164,30 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
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))
Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) ! Compute BSE excitation energies
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)
Bpp(:,:) = Bpp(:,:) + TBs(:,:) - TBt(:,:) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
Cpp(:,:) = Cpp(:,:) + TCs(:,:) - TCt(:,:) allocate(KB_sta(nVVt,nOOt),KC_sta(nVVt,nVVt),KD_sta(nOOt,nOOt))
Dpp(:,:) = Dpp(:,:) + TDs(:,:) - TDt(:,:)
call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin)) if(.not.TDA) 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,eGT,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, &
Taaaa,Tabab,Tbaab,KB_sta)
call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, &
Taaaa,Tabab,Tbaab,KC_sta)
call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, &
Taaaa,Tabab,Tbaab,KD_sta)
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,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)

76
src/GT/RGTpp_ppBSE_static_kernel_B.f90
View File

@ -1,4 +1,4 @@
subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB) subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KB_sta)
! Compute the VVOO block of the static T-matrix ! Compute the VVOO block of the static T-matrix
@ -16,13 +16,11 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: eGF(nBas)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
@ -32,7 +30,11 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
! Output variables ! Output variables
double precision,intent(out) :: TB(nVVx,nOOx) double precision,intent(out) :: KB_sta(nVV,nOO)
! Initialization
KB_sta(:,:) = 0d0
!===============! !===============!
! singlet block ! ! singlet block !
@ -53,16 +55,16 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
chi = 0d0 chi = 0d0
do cd=1,nVV do cd=1,nVV
eps = + Om1(cd) eps = 0d0
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) chi = chi + 0d0
end do end do
do kl=1,nOO do kl=1,nOO
eps = - Om2(kl) eps = 0d0
chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) chi = chi + 0d0
end do end do
TB(ab,ij) = lambda*chi KB_sta(ab,ij) = lambda*chi
end do end do
end do end do
@ -91,54 +93,16 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
chi = 0d0 chi = 0d0
do cd=1,nVV do cd=1,nVV
eps = + Om1(cd) eps = 0d0
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) chi = chi + 0d0
end do end do
do kl=1,nOO do kl=1,nOO
eps = - Om2(kl) eps = 0d0
chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) chi = chi + 0d0
end do end do
TB(ab,ij) = lambda*chi KB_sta(ab,ij) = lambda*chi
end do
end do
end do
end do
end if
!==================!
! alpha-beta block !
!==================!
if(ispin == 3) then
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
ab = ab + 1
ij = 0
do i=nC+1,nO
do j=nC+1,nO
ij = ij + 1
chi = 0d0
do cd=1,nVV
eps = + Om1(cd)
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = - Om2(kl)
chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2)
end do
TB(ab,ij) = lambda*chi
end do end do
end do end do

86
src/GT/RGTpp_ppBSE_static_kernel_C.f90
View File

@ -1,4 +1,4 @@
subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC) subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KC_sta)
! Compute the VVVV block of the static T-matrix ! Compute the VVVV block of the static T-matrix
@ -16,24 +16,27 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: eGF(nBas)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
double precision,external :: Kronecker_delta double precision,external :: Kronecker_delta
double precision :: dem,num
double precision :: chi double precision :: chi
double precision :: eps double precision :: eps
integer :: a,b,c,d,ab,cd,ef,mn integer :: a,b,c,d,ab,cd,ef,mn,m,e
! Output variables ! Output variables
double precision,intent(out) :: TC(nVVx,nVVx) double precision,intent(out) :: KC_sta(nVV,nVV)
! Initialization
KC_sta(:,:) = 0d0
!===============! !===============!
! singlet block ! ! singlet block !
@ -54,18 +57,16 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
chi = 0d0 chi = 0d0
do ef=1,nVV do ef=1,nVV
eps = + Om1(ef) eps = 0d0
chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) & chi = chi + 0d0
+ rho1(a,b,ef)*rho1(d,c,ef)*eps/(eps**2 + eta**2)
end do end do
do mn=1,nOO do mn=1,nOO
eps = - Om2(mn) eps = 0d0
chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) & chi = chi + 0d0
+ rho2(a,b,mn)*rho2(d,c,mn)*eps/(eps**2 + eta**2)
end do end do
TC(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) KC_sta(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
end do end do
end do end do
@ -91,61 +92,18 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
do d=c+1,nBas-nR do d=c+1,nBas-nR
cd = cd + 1 cd = cd + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nBas-nR
do ef=1,nVV dem = eGF(m) - eGF(e)
eps = + Om1(ef) num = 2d0*(Taaaa(a,m,c,e)*Taaaa(e,b,m,d) + Tabab(a,m,c,e)*Tabab(e,b,m,d))
chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2)
end do
do mn=1,nOO
eps = - Om2(mn)
chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2)
end do
TC(ab,cd) = lambda*chi
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
end do end do
end do end do
end do end do
end do end do
end if
!==================!
! alpha-beta block !
!==================!
if(ispin == 3) then
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
cd = cd + 1
chi = 0d0
do ef=1,nVV
eps = + Om1(ef)
chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2)
end do
do mn=1,nOO
eps = - Om2(mn)
chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2)
end do
TC(ab,cd) = lambda*chi
end do
end do
end do end do
end do end do

118
src/GT/RGTpp_ppBSE_static_kernel_D.f90
View File

@ -1,4 +1,4 @@
subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD) subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KD_sta)
! Compute the OOOO block of the static T-matrix ! Compute the OOOO block of the static T-matrix
@ -16,23 +16,25 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
integer,intent(in) :: nOOx
integer,intent(in) :: nVVx
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: eGF(nBas)
double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: rho2(nBas,nBas,nOO) double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
! Local variables
double precision,external :: Kronecker_delta
double precision :: dem,num
double precision :: chi double precision :: chi
double precision :: eps double precision :: eps
integer :: i,j,k,l,ij,kl,ef,mn integer :: i,j,k,l,ij,kl,ef,mn,m,e
! Output variables ! Output variables
double precision,intent(out) :: TD(nOOx,nOOx) double precision,intent(out) :: KD_sta(nOO,nOO)
! Initialization
KD_sta(:,:) = 0d0
!===============! !===============!
! singlet block ! ! singlet block !
@ -50,20 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
do l=k,nO do l=k,nO
kl = kl + 1 kl = kl + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nBas-nR
dem = eGF(m) - eGF(e)
! Wabab_{ijkl}
num = Taaaa(i,m,k,e)*Tabab(e,j,m,l) + Tabab(i,m,k,e)*Taaaa(e,j,m,l)
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
do ef=1,nVV num = Taaaa(i,e,k,m)*Tabab(m,j,e,l) + Tabab(i,e,k,m)*Taaaa(m,j,e,l)
eps = + Om1(ef) KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
num = Taaaa(j,m,k,e)*Tabab(e,i,m,l) + Tabab(j,m,k,e)*Taaaa(e,i,m,l)
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
num = Taaaa(j,e,k,m)*Tabab(m,i,e,l) + Tabab(j,e,k,m)*Taaaa(m,i,e,l)
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
num = Tabab(i,m,k,e)*Tbaab(e,j,m,l) + Tbaab(i,e,k,m)*Tabab(m,j,e,l)
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
num = Tabab(j,m,k,e)*Tbaab(e,i,m,l) + Tbaab(j,e,k,m)*Tabab(m,i,e,l)
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
end do
end do end do
do mn=1,nOO KD_sta(ij,kl) = KD_sta(ij,kl) / sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
eps = - Om2(mn)
chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2)
end do
TD(ij,kl) = lambda*chi
end do end do
end do end do
@ -88,19 +101,24 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
do l=k+1,nO do l=k+1,nO
kl = kl + 1 kl = kl + 1
chi = 0d0 do m=nC+1,nO
do e=nO+1,nBas-nR
dem = eGF(m) - eGF(e)
num = Taaaa(i,m,k,e)*Taaaa(e,j,m,l) + Tabab(i,m,k,e)*Tabab(e,j,m,l)
do ef=1,nVV KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
eps = + Om1(ef)
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
end do
do mn=1,nOO num = Taaaa(i,e,k,m)*Taaaa(m,j,e,l) + Tabab(i,e,k,m)*Tabab(m,j,e,l)
eps = - Om2(mn)
chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2)
end do
TD(ij,kl) = lambda*chi KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
num = Taaaa(j,m,k,e)*Taaaa(e,i,m,l) + Tabab(j,m,k,e)*Tabab(e,i,m,l)
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
num = Taaaa(j,e,k,m)*Taaaa(m,i,e,l) + Tabab(j,e,k,m)*Tabab(m,i,e,l)
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
end do end do
end do end do
@ -108,45 +126,9 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
end do end do
end do end do
end if
!==================!
! alpha-beta block !
!==================!
if(ispin == 3) then
ij = 0
do i=nC+1,nO
do j=nC+1,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=nC+1,nO
kl = kl + 1
chi = 0d0
do ef=1,nVV
eps = + Om1(ef)
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
end do
do mn=1,nOO
eps = - Om2(mn)
chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2)
end do
TD(ij,kl) = lambda*chi
end do end do
end do end do
end do
end do
end if end if
end subroutine end subroutine

16
src/GT/RGTpp_self_energy.f90
View File

@ -48,14 +48,14 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho
do cd=1,nVVs do cd=1,nVVs
eps = e(p) + e(i) - Om1s(cd) eps = e(p) + e(i) - Om1s(cd)
num = rho1s(p,i,cd)*rho1s(q,i,cd) num = (1d0/2d0)*rho1s(p,i,cd)*rho1s(q,i,cd)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = e(p) + e(i) - Om1t(cd) eps = e(p) + e(i) - Om1t(cd)
num = rho1t(p,i,cd)*rho1t(q,i,cd) num = (3d0/2d0)*rho1t(p,i,cd)*rho1t(q,i,cd)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
@ -74,14 +74,14 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho
do kl=1,nOOs do kl=1,nOOs
eps = e(p) + e(a) - Om2s(kl) eps = e(p) + e(a) - Om2s(kl)
num = rho2s(p,a,kl)*rho2s(q,a,kl) num = (1d0/2d0)*rho2s(p,a,kl)*rho2s(q,a,kl)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = e(p) + e(a) - Om2t(kl) eps = e(p) + e(a) - Om2t(kl)
num = rho2t(p,a,kl)*rho2t(q,a,kl) num = (3d0/2d0)*rho2t(p,a,kl)*rho2t(q,a,kl)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
@ -99,13 +99,13 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho
do cd=1,nVVs do cd=1,nVVs
eps = e(i) + e(j) - Om1s(cd) eps = e(i) + e(j) - Om1s(cd)
num = rho1s(i,j,cd)*rho1s(i,j,cd) num = (1d0/2d0)*rho1s(i,j,cd)*rho1s(i,j,cd)
EcGM = EcGM + num*eps/(eps**2 + eta**2) EcGM = EcGM + num*eps/(eps**2 + eta**2)
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = e(i) + e(j) - Om1t(cd) eps = e(i) + e(j) - Om1t(cd)
num = rho1t(i,j,cd)*rho1t(i,j,cd) num = (3d0/2d0)*rho1t(i,j,cd)*rho1t(i,j,cd)
EcGM = EcGM + num*eps/(eps**2 + eta**2) EcGM = EcGM + num*eps/(eps**2 + eta**2)
end do end do
@ -117,13 +117,13 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho
do kl=1,nOOs do kl=1,nOOs
eps = e(a) + e(b) - Om2s(kl) eps = e(a) + e(b) - Om2s(kl)
num = rho2s(a,b,kl)*rho2s(a,b,kl) num = (1d0/2d0)*rho2s(a,b,kl)*rho2s(a,b,kl)
EcGM = EcGM - num*eps/(eps**2 + eta**2) EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = e(a) + e(b) - Om2t(kl) eps = e(a) + e(b) - Om2t(kl)
num = rho2t(a,b,kl)*rho2t(a,b,kl) num = (3d0/2d0)*rho2t(a,b,kl)*rho2t(a,b,kl)
EcGM = EcGM - num*eps/(eps**2 + eta**2) EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do

16
src/GT/RGTpp_self_energy_diag.f90
View File

@ -48,14 +48,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1
do cd=1,nVVs do cd=1,nVVs
eps = e(p) + e(i) - Om1s(cd) eps = e(p) + e(i) - Om1s(cd)
num = rho1s(p,i,cd)**2 num = (1d0/2d0)*rho1s(p,i,cd)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = e(p) + e(i) - Om1t(cd) eps = e(p) + e(i) - Om1t(cd)
num = rho1t(p,i,cd)**2 num = (3d0/2d0)*rho1t(p,i,cd)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
@ -72,14 +72,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1
do kl=1,nOOs do kl=1,nOOs
eps = e(p) + e(a) - Om2s(kl) eps = e(p) + e(a) - Om2s(kl)
num = rho2s(p,a,kl)**2 num = (1d0/2d0)*rho2s(p,a,kl)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = e(p) + e(a) - Om2t(kl) eps = e(p) + e(a) - Om2t(kl)
num = rho2t(p,a,kl)**2 num = (3d0/2d0)*rho2t(p,a,kl)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
@ -96,13 +96,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1
do cd=1,nVVs do cd=1,nVVs
eps = e(i) + e(j) - Om1s(cd) eps = e(i) + e(j) - Om1s(cd)
num = rho1s(i,j,cd)**2 num = (1d0/2d0)*rho1s(i,j,cd)**2
EcGM = EcGM + num*eps/(eps**2 + eta**2) EcGM = EcGM + num*eps/(eps**2 + eta**2)
end do end do
do cd=1,nVVt do cd=1,nVVt
eps = e(i) + e(j) - Om1t(cd) eps = e(i) + e(j) - Om1t(cd)
num = rho1t(i,j,cd)**2 num = (3d0/2d0)*rho1t(i,j,cd)**2
EcGM = EcGM + num*eps/(eps**2 + eta**2) EcGM = EcGM + num*eps/(eps**2 + eta**2)
end do end do
@ -114,13 +114,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1
do kl=1,nOOs do kl=1,nOOs
eps = e(a) + e(b) - Om2s(kl) eps = e(a) + e(b) - Om2s(kl)
num = rho2s(a,b,kl)**2 num = (1d0/2d0)*rho2s(a,b,kl)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2) EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do
do kl=1,nOOt do kl=1,nOOt
eps = e(a) + e(b) - Om2t(kl) eps = e(a) + e(b) - Om2t(kl)
num = rho2t(a,b,kl)**2 num = (3d0/2d0)*rho2t(a,b,kl)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2) EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do end do

2
src/GW/evRGW.f90
View File

@ -85,7 +85,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! SRG regularization ! SRG regularization
flow = 500d0 flow = 100d0
if(doSRG) then if(doSRG) then

59
src/LR/ppGLR.f90
View File

@ -49,7 +49,6 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP)) allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP))
! Hermitian case for TDA ! Hermitian case for TDA
if(TDA) then if(TDA) then
X1(:,:) = +Cpp(:,:) X1(:,:) = +Cpp(:,:)
@ -72,7 +71,7 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
M( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO) M( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO)
M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO)) M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO))
! if((nOO .eq. 0) .or. (nVV .eq. 0)) then ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! Diagonalize the pp matrix ! Diagonalize the pp matrix
@ -82,36 +81,36 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! else ! else
! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 ! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 ! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not ! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
! imp_bio = .True. ! impose bi-orthogonality ! imp_bio = .True. ! impose bi-orthogonality
! verbose = .False. ! verbose = .False.
! call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) ! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
!
! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
!
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! endif ! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! endif
end if end if

56
src/LR/ppLR.f90
View File

@ -82,36 +82,36 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! else ! else
! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 ! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 ! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not ! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
! imp_bio = .True. ! impose bi-orthogonality ! imp_bio = .True. ! impose bi-orthogonality
! verbose = .False. ! verbose = .False.
! call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) ! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
!
! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
!
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! endif ! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
! endif
end if end if

33
src/QuAcK/GQuAcK.f90
View File

@ -1,11 +1,12 @@
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
maxSCF_CC,max_diis_CC,thresh_CC, & maxSCF_CC,max_diis_CC,thresh_CC, &
TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
implicit none implicit none
@ -23,6 +24,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA
logical,intent(in) :: doG0F2,doevGF2,doqsGF2 logical,intent(in) :: doG0F2,doevGF2,doqsGF2
logical,intent(in) :: doG0W0,doevGW,doqsGW logical,intent(in) :: doG0W0,doevGW,doqsGW
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp
integer,intent(in) :: nNuc,nBas integer,intent(in) :: nNuc,nBas
integer,intent(in) :: nC integer,intent(in) :: nC
@ -62,12 +64,17 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
logical,intent(in) :: TDA_W,lin_GW,reg_GW logical,intent(in) :: TDA_W,lin_GW,reg_GW
double precision,intent(in) :: eta_GW double precision,intent(in) :: eta_GW
integer,intent(in) :: maxSCF_GT,max_diis_GT
double precision,intent(in) :: thresh_GT
logical,intent(in) :: TDA_T,lin_GT,reg_GT
double precision,intent(in) :: eta_GT
logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA
logical,intent(in) :: doACFDT,exchange_kernel,doXBS logical,intent(in) :: doACFDT,exchange_kernel,doXBS
! Local variables ! Local variables
logical :: doMP,doCC,doRPA,doGF,doGW logical :: doMP,doCC,doRPA,doGF,doGW,doGT
double precision :: start_HF ,end_HF ,t_HF double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab double precision :: start_stab ,end_stab ,t_stab
@ -77,6 +84,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_GF ,end_GF ,t_GF double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GW ,end_GW ,t_GW double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT
double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:) double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:)
double precision :: EGHF double precision :: EGHF
@ -290,4 +298,25 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
end if end if
!-----------------!
! T-matrix module !
!-----------------!
doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp
if(doGT) then
call wall_time(start_GT)
call GGT(dotest,doG0T0pp,doevGTpp,doqsGTpp, &
maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
TDA_T,TDA,dBSE,dTDA,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO, &
dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds'
write(*,*)
end if
end subroutine end subroutine

4
src/QuAcK/QuAcK.f90
View File

@ -237,15 +237,15 @@ program QuAcK
!--------------------------! !--------------------------!
! Generalized QuAcK branch ! ! Generalized QuAcK branch !
!--------------------------! !--------------------------!
if(doGQuAcK) & if(doGQuAcK) &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, &
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
!-----------! !-----------!