mirror of
https://github.com/pfloos/quack
synced 2025-01-03 01:55:57 +01:00
Merge branch 'master' into work_env
This commit is contained in:
Abdallah Ammar
commit
e4f89c21f5
@ -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)
|
||||
|
||||
! 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) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: dophBSE,dophBSE2
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA_T
|
||||
logical,intent(in) :: TDA
|
||||
@ -37,7 +37,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: print_T = .false.
|
||||
logical :: print_T = .true.
|
||||
|
||||
integer :: nOO
|
||||
integer :: nVV
|
||||
@ -72,7 +72,6 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
|
||||
write(*,*)'**********************************'
|
||||
write(*,*)
|
||||
|
||||
|
||||
! TDA for T
|
||||
|
||||
if(TDA_T) then
|
||||
@ -104,7 +103,9 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
|
||||
! Compute linear response
|
||||
|
||||
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_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)
|
||||
@ -114,12 +115,11 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
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
|
||||
!----------------------------------------------
|
||||
|
||||
call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
|
||||
!----------------------------------------------
|
||||
@ -170,6 +170,23 @@ 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)
|
||||
|
||||
!----------------------------------------------
|
||||
! 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@GHF correlation energy = ',EcBSE,' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
117
src/GT/GGT.f90
Normal file
117
src/GT/GGT.f90
Normal 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
|
||||
66
src/GT/GGT_Tmatrix.f90
Normal file
66
src/GT/GGT_Tmatrix.f90
Normal file
@ -0,0 +1,66 @@
|
||||
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 :: 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,nOO,nVV,T,ERI,rho1,rho2,Om1,Om2) &
|
||||
!$OMP PRIVATE(p,q,r,s,cd,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)
|
||||
|
||||
do cd=1,nVV
|
||||
T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd)*rho1(r,s,cd)/Om1(cd)
|
||||
end do
|
||||
|
||||
do kl=1,nOO
|
||||
T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl)*rho2(r,s,kl)/Om2(kl)
|
||||
end do
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
end subroutine
|
||||
@ -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
|
||||
|
||||
@ -12,7 +12,6 @@ subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
|
||||
124
src/GT/GGTpp_ppBSE.f90
Normal file
124
src/GT/GGTpp_ppBSE.f90
Normal file
@ -0,0 +1,124 @@
|
||||
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
|
||||
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))
|
||||
|
||||
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
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)
|
||||
|
||||
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))
|
||||
|
||||
if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
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)
|
||||
|
||||
!----------------------------------------------
|
||||
! 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))
|
||||
|
||||
if(.not.TDA) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta)
|
||||
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)
|
||||
|
||||
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,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
|
||||
@ -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
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
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(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
|
||||
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
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,a,b,ij,ab,cd,kl
|
||||
double precision :: dem,num
|
||||
integer :: p,q,r,s,e,m
|
||||
integer :: i,j,a,b
|
||||
integer :: ij,ab,kl,cd
|
||||
|
||||
! 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
|
||||
do a=nO+1,nBas-nR
|
||||
do b=a+1,nBas-nR
|
||||
do a=nO+1,nOrb-nR
|
||||
do b=a+1,nOrb-nR
|
||||
ab = ab + 1
|
||||
|
||||
ij = 0
|
||||
@ -44,19 +50,13 @@ subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nOrb-nR
|
||||
|
||||
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
|
||||
dem = eGF(m) - eGF(e)
|
||||
num = T(a,m,i,e) * T(e,b,m,j) + T(a,e,i,m) * T(m,b,e,j)
|
||||
num = num - T(b,m,i,e) * T(e,a,m,j) - T(b,e,i,m) * T(m,a,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -64,5 +64,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 subroutine
|
||||
|
||||
@ -1,64 +1,67 @@
|
||||
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
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
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(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
|
||||
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
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: a,b,c,d,ab,cd,ef,mn
|
||||
double precision :: dem,num
|
||||
integer :: p,q,r,s,e,m
|
||||
integer :: a,b,c,d,k,l
|
||||
integer :: ab,kl,cd
|
||||
|
||||
! 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
|
||||
do a=nO+1,nBas-nR
|
||||
do b=a+1,nBas-nR
|
||||
do a=nO+1,nOrb-nR
|
||||
do b=a+1,nOrb-nR
|
||||
ab = ab + 1
|
||||
|
||||
cd = 0
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c+1,nBas-nR
|
||||
do c=nO+1,nOrb-nR
|
||||
do d=c+1,nOrb-nR
|
||||
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(e,b,m,d) + T(a,e,c,m) * T(m,b,e,d)
|
||||
num = num - T(b,m,c,e) * T(e,a,m,d) - T(b,e,c,m) * T(m,a,e,d)
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
|
||||
@ -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
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
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(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: rho1(nOrb,nOrb,nVV)
|
||||
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
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,k,l,ij,kl,ef,mn
|
||||
double precision :: dem,num
|
||||
integer :: p,q,r,s,e,m
|
||||
integer :: i,j,k,l
|
||||
integer :: ij,kl,cd
|
||||
|
||||
! 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
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
@ -44,24 +50,21 @@ subroutine GGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nOrb-nR
|
||||
|
||||
do ef=1,nVV
|
||||
eps = + Om1(ef)
|
||||
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
dem = eGF(m) - eGF(e)
|
||||
num = T(i,m,k,e) * T(e,j,m,l) + T(i,e,k,m) * T(m,j,e,l)
|
||||
num = num - T(j,m,k,e) * T(e,i,m,l) - T(j,e,k,m) * T(m,i,e,l)
|
||||
|
||||
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 subroutine
|
||||
|
||||
@ -81,7 +81,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh
|
||||
|
||||
do cd=1,nVV
|
||||
eps = e(i) + e(j) - Om1(cd)
|
||||
num = rho1(i,j,cd)**2
|
||||
num = 0.5d0*rho1(i,j,cd)**2
|
||||
EcGM = EcGM + num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
@ -93,7 +93,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh
|
||||
|
||||
do kl=1,nOO
|
||||
eps = e(a) + e(b) - Om2(kl)
|
||||
num = rho2(a,b,kl)**2
|
||||
num = 0.5d0*rho2(a,b,kl)**2
|
||||
EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
|
||||
@ -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(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)
|
||||
|
||||
@ -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(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)
|
||||
|
||||
|
||||
@ -1,6 +1,3 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
|
||||
maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
|
||||
doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, &
|
||||
|
||||
157
src/GT/RGT_Tmatrix.f90
Normal file
157
src/GT/RGT_Tmatrix.f90
Normal file
@ -0,0 +1,157 @@
|
||||
subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,T)
|
||||
|
||||
! Compute the T-matrix tensor elements
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOOs, nOOt
|
||||
integer,intent(in) :: nVVs, nVVt
|
||||
integer,intent(in) :: isp_T
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Om1s(nVVs)
|
||||
double precision,intent(in) :: rho1s(nBas,nBas,nVVs)
|
||||
double precision,intent(in) :: Om2s(nOOs)
|
||||
double precision,intent(in) :: rho2s(nBas,nBas,nOOs)
|
||||
double precision,intent(in) :: Om1t(nVVt)
|
||||
double precision,intent(in) :: rho1t(nBas,nBas,nVVt)
|
||||
double precision,intent(in) :: Om2t(nOOt)
|
||||
double precision,intent(in) :: rho2t(nBas,nBas,nOOt)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
integer :: p,q,r,s
|
||||
integer :: c,d,k,l
|
||||
integer :: kl,cd
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: T(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Initialization
|
||||
|
||||
T(:,:,:,:) = 0d0
|
||||
|
||||
! Elements aaaa
|
||||
|
||||
if(isp_T == 1) then
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOt,nVVt,rho1t,rho2t,Om1t,Om2t) &
|
||||
!$OMP PRIVATE(p,q,r,s,cd,kl) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do s=nC+1,nBas-nR
|
||||
do r=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r)
|
||||
|
||||
do cd=1,nVVt
|
||||
T(p,q,r,s) = T(p,q,r,s) - rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd)
|
||||
end do ! cd
|
||||
|
||||
do kl=1,nOOt
|
||||
T(p,q,r,s) = T(p,q,r,s) + rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl)
|
||||
enddo ! kl
|
||||
|
||||
enddo ! p
|
||||
enddo ! q
|
||||
enddo ! r
|
||||
enddo ! s
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
! Elements abab
|
||||
|
||||
if(isp_T == 2) then
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOs,nOOt,nVVs,nVVt,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) &
|
||||
!$OMP PRIVATE(p,q,r,s,cd,kl) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do s=nC+1,nBas-nR
|
||||
do r=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
T(p,q,r,s) = ERI(p,q,r,s)
|
||||
|
||||
do cd=1,nVVs
|
||||
T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1s(p,q,cd) * rho1s(r,s,cd) / Om1s(cd)
|
||||
end do ! cd
|
||||
|
||||
do cd=1,nVVt
|
||||
T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd)
|
||||
end do ! cd
|
||||
|
||||
do kl=1,nOOs
|
||||
T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2s(p,q,kl) * rho2s(r,s,kl) / Om2s(kl)
|
||||
enddo ! kl
|
||||
|
||||
do kl=1,nOOt
|
||||
T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl)
|
||||
enddo ! kl
|
||||
|
||||
enddo ! p
|
||||
enddo ! q
|
||||
enddo ! r
|
||||
enddo ! s
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
! Elements baab
|
||||
|
||||
if(isp_T == 3) then
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOs,nOOt,nVVs,nVVt,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) &
|
||||
!$OMP PRIVATE(p,q,r,s,cd,kl) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do s=nC+1,nBas-nR
|
||||
do r=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
T(p,q,r,s) = - ERI(p,q,s,r)
|
||||
|
||||
do cd=1,nVVs
|
||||
T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho1s(p,q,cd) * rho1s(r,s,cd) / Om1s(cd)
|
||||
end do ! cd
|
||||
|
||||
do cd=1,nVVt
|
||||
T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd)
|
||||
end do ! cd
|
||||
|
||||
do kl=1,nOOs
|
||||
T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho2s(p,q,kl) * rho2s(r,s,kl) / Om2s(kl)
|
||||
enddo ! kl
|
||||
|
||||
do kl=1,nOOt
|
||||
T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl)
|
||||
enddo ! kl
|
||||
|
||||
enddo ! p
|
||||
enddo ! q
|
||||
enddo ! r
|
||||
enddo ! s
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
@ -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
|
||||
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
|
||||
|
||||
do cd=1,nVVt
|
||||
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
|
||||
@ -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
|
||||
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
|
||||
|
||||
do kl=1,nOOt
|
||||
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
|
||||
|
||||
@ -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
|
||||
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
|
||||
|
||||
do cd=1,nVVt
|
||||
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
|
||||
|
||||
!----------------------------------------------
|
||||
! Virtual part of the T-matrix self-energy
|
||||
!
|
||||
!----------------------------------------------
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
do kl=1,nOOs
|
||||
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
|
||||
|
||||
do kl=1,nOOt
|
||||
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
|
||||
|
||||
@ -66,7 +66,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
cd = cd + 1
|
||||
rho1(p,q,ab) = rho1(p,q,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
|
||||
|
||||
@ -76,7 +76,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
kl = kl + 1
|
||||
rho1(p,q,ab) = rho1(p,q,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
|
||||
|
||||
@ -93,7 +93,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
cd = cd + 1
|
||||
rho2(p,q,ij) = rho2(p,q,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
|
||||
|
||||
@ -103,7 +103,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
kl = kl + 1
|
||||
rho2(p,q,ij) = rho2(p,q,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
|
||||
|
||||
@ -148,7 +148,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
cd = cd + 1
|
||||
|
||||
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 ! c
|
||||
|
||||
@ -159,7 +159,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
kl = kl + 1
|
||||
|
||||
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 ! k
|
||||
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
|
||||
|
||||
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 ! c
|
||||
|
||||
@ -190,7 +190,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
|
||||
kl = kl + 1
|
||||
|
||||
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 ! k
|
||||
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)
|
||||
|
||||
rho1 = rho1*sqrt(3d0/2d0)
|
||||
rho2 = rho2*sqrt(3d0/2d0)
|
||||
rho1 = rho1
|
||||
rho2 = rho2
|
||||
|
||||
endif
|
||||
endif
|
||||
|
||||
@ -34,7 +34,7 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
integer :: ispin, isp_T
|
||||
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
|
||||
@ -46,66 +46,89 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
|
||||
double precision,allocatable :: X2s(:,:),X2t(:,:)
|
||||
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
|
||||
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
|
||||
double precision,allocatable :: TBs(:,:),TCs(:,:),TDs(:,:)
|
||||
double precision,allocatable :: TBt(:,:),TCt(:,:),TDt(:,:)
|
||||
double precision,allocatable :: KB_sta(:,:)
|
||||
double precision,allocatable :: KC_sta(:,:)
|
||||
double precision,allocatable :: KD_sta(:,:)
|
||||
double precision,allocatable :: Taaaa(:,:,:,:),Tabab(:,:,:,:),Tbaab(:,:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcBSE(nspin)
|
||||
|
||||
!------------------------------------!
|
||||
! Compute T-matrix for singlet block !
|
||||
!------------------------------------!
|
||||
!---------------------------------
|
||||
! Compute ppRPA excitation density
|
||||
!---------------------------------
|
||||
|
||||
ispin = 1
|
||||
! Singlet contribution
|
||||
|
||||
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))
|
||||
isp_T = 1
|
||||
|
||||
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)
|
||||
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
|
||||
|
||||
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
|
||||
if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
|
||||
call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
allocate(TBs(nVVs,nOOs),TCs(nVVs,nVVs),TDs(nOOs,nOOs))
|
||||
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs))
|
||||
allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(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)
|
||||
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(Om1s,X1s,Y1s,Om2s,X2s,Y2s)
|
||||
allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs))
|
||||
|
||||
!------------------------------------!
|
||||
! Compute T-matrix for triplet block !
|
||||
!------------------------------------!
|
||||
call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
ispin = 2
|
||||
deallocate(X1s,Y1s,X2s,Y2s,Bpp,Cpp,Dpp)
|
||||
|
||||
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))
|
||||
! Triplet contribution
|
||||
|
||||
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)
|
||||
isp_T = 2
|
||||
|
||||
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
|
||||
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
allocate(TBt(nVVt,nOOt),TCt(nVVt,nVVt),TDt(nOOt,nOOt))
|
||||
if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
|
||||
call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp)
|
||||
call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp)
|
||||
|
||||
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)
|
||||
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt))
|
||||
allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
|
||||
|
||||
deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t)
|
||||
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(rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
|
||||
|
||||
call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
deallocate(X1t,Y1t,X2t,Y2t,Bpp,Cpp,Dpp)
|
||||
|
||||
!---------------------------------
|
||||
! Compute T matrix elements
|
||||
!---------------------------------
|
||||
|
||||
! Elements aaaa
|
||||
|
||||
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)
|
||||
|
||||
! Elements abab
|
||||
|
||||
isp_T = 2
|
||||
allocate(Tabab(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,Tabab)
|
||||
|
||||
! Elements baab
|
||||
|
||||
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 !
|
||||
@ -115,22 +138,31 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
|
||||
|
||||
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))
|
||||
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs))
|
||||
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)
|
||||
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(:,:)
|
||||
Cpp(:,:) = Cpp(:,:) - TCs(:,:) - TCt(:,:)
|
||||
Dpp(:,:) = Dpp(:,:) - TDs(:,:) - TDt(:,:)
|
||||
if(.not.TDA) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KB_sta)
|
||||
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_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
|
||||
|
||||
@ -144,18 +176,27 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
|
||||
|
||||
EcBSE(ispin) = 0d0
|
||||
|
||||
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))
|
||||
allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt))
|
||||
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)
|
||||
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)
|
||||
! Compute BSE excitation energies
|
||||
|
||||
Bpp(:,:) = Bpp(:,:) + TBs(:,:) - TBt(:,:)
|
||||
Cpp(:,:) = Cpp(:,:) + TCs(:,:) - TCt(:,:)
|
||||
Dpp(:,:) = Dpp(:,:) + TDs(:,:) - TDt(:,:)
|
||||
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
|
||||
allocate(KB_sta(nVVt,nOOt),KC_sta(nVVt,nVVt),KD_sta(nOOt,nOOt))
|
||||
|
||||
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) 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)
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -16,23 +16,25 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
double precision,intent(in) :: eGF(nBas)
|
||||
double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,a,b,ij,ab,cd,kl
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: dem,num
|
||||
integer :: i,j,a,b,ij,ab,cd,kl,m,e
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: TB(nVVx,nOOx)
|
||||
double precision,intent(out) :: KB_sta(nVV,nOO)
|
||||
|
||||
! Initialization
|
||||
|
||||
KB_sta(:,:) = 0d0
|
||||
|
||||
!===============!
|
||||
! singlet block !
|
||||
@ -50,19 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do j=i,nO
|
||||
ij = ij + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nBas-nR
|
||||
dem = eGF(m) - eGF(e)
|
||||
! Wabab_{ijkl}
|
||||
num = Taaaa(a,m,i,e)*Tabab(e,b,m,j) + Tabab(a,m,i,e)*Taaaa(e,b,m,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
do cd=1,nVV
|
||||
eps = + Om1(cd)
|
||||
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
|
||||
num = Taaaa(a,e,i,m)*Tabab(m,b,e,j) + Tabab(a,e,i,m)*Taaaa(m,b,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,m,i,e)*Tabab(e,a,m,j) + Tabab(b,m,i,e)*Taaaa(e,a,m,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,e,i,m)*Tabab(m,a,e,j) + Tabab(b,e,i,m)*Taaaa(m,a,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Tbaab(a,m,i,e)*Tbaab(e,b,m,j) + Tbaab(a,e,i,m)*Tbaab(m,b,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Tbaab(b,m,i,e)*Tbaab(e,a,m,j) + Tbaab(b,e,i,m)*Tbaab(m,a,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
end do
|
||||
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
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) / sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -88,19 +102,21 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nBas-nR
|
||||
dem = eGF(m) - eGF(e)
|
||||
|
||||
do cd=1,nVV
|
||||
eps = + Om1(cd)
|
||||
chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(a,m,i,e)*Taaaa(e,b,m,j) + Tabab(a,m,i,e)*Tabab(e,b,m,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
do kl=1,nOO
|
||||
eps = - Om2(kl)
|
||||
chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(a,e,i,m)*Taaaa(m,b,e,j) + Tabab(a,e,i,m)*Tabab(m,b,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
TB(ab,ij) = lambda*chi
|
||||
num = Taaaa(b,m,i,e)*Taaaa(e,a,m,j) + Tabab(b,m,i,e)*Tabab(e,a,m,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,e,i,m)*Taaaa(m,a,e,j) + Tabab(b,e,i,m)*Tabab(m,a,e,j)
|
||||
KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -108,41 +124,6 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -16,24 +16,25 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
double precision,intent(in) :: eGF(nBas)
|
||||
double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: a,b,c,d,ab,cd,ef,mn
|
||||
double precision :: dem,num
|
||||
integer :: a,b,c,d,ab,cd,ef,mn,m,e
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: TC(nVVx,nVVx)
|
||||
double precision,intent(out) :: KC_sta(nVV,nVV)
|
||||
|
||||
! Initialization
|
||||
|
||||
KC_sta(:,:) = 0d0
|
||||
|
||||
!===============!
|
||||
! singlet block !
|
||||
@ -51,21 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do d=c,nBas-nR
|
||||
cd = cd + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nBas-nR
|
||||
dem = eGF(m) - eGF(e)
|
||||
! Wabab_{ijkl}
|
||||
num = Taaaa(a,m,c,e)*Tabab(e,b,m,d) + Tabab(a,m,c,e)*Taaaa(e,b,m,d)
|
||||
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) &
|
||||
+ rho1(a,b,ef)*rho1(d,c,ef)*eps/(eps**2 + eta**2)
|
||||
num = Taaaa(a,e,c,m)*Tabab(m,b,e,d) + Tabab(a,e,c,m)*Taaaa(m,b,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,m,c,e)*Tabab(e,a,m,d) + Tabab(b,m,c,e)*Taaaa(e,a,m,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,e,c,m)*Tabab(m,a,e,d) + Tabab(b,e,c,m)*Taaaa(m,a,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Tbaab(a,m,c,e)*Tbaab(e,b,m,d) + Tbaab(a,e,c,m)*Tbaab(m,b,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Tbaab(b,m,c,e)*Tbaab(e,a,m,d) + Tbaab(b,e,c,m)*Tbaab(m,a,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**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) &
|
||||
+ rho2(a,b,mn)*rho2(d,c,mn)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
TC(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) / sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -91,19 +102,21 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do d=c+1,nBas-nR
|
||||
cd = cd + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nBas-nR
|
||||
dem = eGF(m) - eGF(e)
|
||||
|
||||
do ef=1,nVV
|
||||
eps = + Om1(ef)
|
||||
chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(a,m,c,e)*Taaaa(e,b,m,d) + Tabab(a,m,c,e)*Tabab(e,b,m,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
do mn=1,nOO
|
||||
eps = - Om2(mn)
|
||||
chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(a,e,c,m)*Taaaa(m,b,e,d) + Tabab(a,e,c,m)*Tabab(m,b,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
TC(ab,cd) = lambda*chi
|
||||
num = Taaaa(b,m,c,e)*Taaaa(e,a,m,d) + Tabab(b,m,c,e)*Tabab(e,a,m,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Taaaa(b,e,c,m)*Taaaa(m,a,e,d) + Tabab(b,e,c,m)*Tabab(m,a,e,d)
|
||||
KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -111,41 +124,6 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -16,23 +16,23 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
integer,intent(in) :: nOOx
|
||||
integer,intent(in) :: nVVx
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
double precision,intent(in) :: eGF(nBas)
|
||||
double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,k,l,ij,kl,ef,mn
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: dem,num
|
||||
integer :: i,j,k,l,ij,kl,ef,mn,m,e
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: TD(nOOx,nOOx)
|
||||
double precision,intent(out) :: KD_sta(nOO,nOO)
|
||||
|
||||
! Initialization
|
||||
|
||||
KD_sta(:,:) = 0d0
|
||||
|
||||
!===============!
|
||||
! singlet block !
|
||||
@ -50,20 +50,31 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do l=k,nO
|
||||
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
|
||||
eps = + Om1(ef)
|
||||
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
|
||||
num = Taaaa(i,e,k,m)*Tabab(m,j,e,l) + Tabab(i,e,k,m)*Taaaa(m,j,e,l)
|
||||
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**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 = Tbaab(i,m,k,e)*Tbaab(e,j,m,l) + Tbaab(i,e,k,m)*Tbaab(m,j,e,l)
|
||||
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
num = Tbaab(j,m,k,e)*Tbaab(e,i,m,l) + Tbaab(j,e,k,m)*Tbaab(m,i,e,l)
|
||||
KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
|
||||
end do
|
||||
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) / sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -88,19 +99,21 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
chi = 0d0
|
||||
do m=nC+1,nO
|
||||
do e=nO+1,nBas-nR
|
||||
dem = eGF(m) - eGF(e)
|
||||
|
||||
do ef=1,nVV
|
||||
eps = + Om1(ef)
|
||||
chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(i,m,k,e)*Taaaa(e,j,m,l) + Tabab(i,m,k,e)*Tabab(e,j,m,l)
|
||||
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
do mn=1,nOO
|
||||
eps = - Om2(mn)
|
||||
chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
num = Taaaa(i,e,k,m)*Taaaa(m,j,e,l) + Tabab(i,e,k,m)*Tabab(m,j,e,l)
|
||||
KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
TD(ij,kl) = lambda*chi
|
||||
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
|
||||
@ -108,45 +121,9 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n
|
||||
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 if
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -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
|
||||
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)
|
||||
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
do cd=1,nVVt
|
||||
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)
|
||||
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
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
|
||||
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)
|
||||
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
do kl=1,nOOt
|
||||
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)
|
||||
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
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
|
||||
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)
|
||||
end do
|
||||
|
||||
do cd=1,nVVt
|
||||
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)
|
||||
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
|
||||
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)
|
||||
end do
|
||||
|
||||
do kl=1,nOOt
|
||||
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)
|
||||
end do
|
||||
|
||||
|
||||
@ -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
|
||||
eps = e(p) + e(i) - Om1s(cd)
|
||||
num = rho1s(p,i,cd)**2
|
||||
num = 0.5d0*rho1s(p,i,cd)**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
|
||||
end do
|
||||
|
||||
do cd=1,nVVt
|
||||
eps = e(p) + e(i) - Om1t(cd)
|
||||
num = rho1t(p,i,cd)**2
|
||||
num = 1.5d0*rho1t(p,i,cd)**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
|
||||
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
|
||||
eps = e(p) + e(a) - Om2s(kl)
|
||||
num = rho2s(p,a,kl)**2
|
||||
num = 0.5d0*rho2s(p,a,kl)**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
|
||||
end do
|
||||
|
||||
do kl=1,nOOt
|
||||
eps = e(p) + e(a) - Om2t(kl)
|
||||
num = rho2t(p,a,kl)**2
|
||||
num = 1.5d0*rho2t(p,a,kl)**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
|
||||
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
|
||||
eps = e(i) + e(j) - Om1s(cd)
|
||||
num = rho1s(i,j,cd)**2
|
||||
num = 0.5d0*rho1s(i,j,cd)**2
|
||||
EcGM = EcGM + num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
do cd=1,nVVt
|
||||
eps = e(i) + e(j) - Om1t(cd)
|
||||
num = rho1t(i,j,cd)**2
|
||||
num = 1.5d0*rho1t(i,j,cd)**2
|
||||
EcGM = EcGM + num*eps/(eps**2 + eta**2)
|
||||
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
|
||||
eps = e(a) + e(b) - Om2s(kl)
|
||||
num = rho2s(a,b,kl)**2
|
||||
num = 0.5d0*rho2s(a,b,kl)**2
|
||||
EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
do kl=1,nOOt
|
||||
eps = e(a) + e(b) - Om2t(kl)
|
||||
num = rho2t(a,b,kl)**2
|
||||
num = 1.5d0*rho2t(a,b,kl)**2
|
||||
EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
|
||||
subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,EGHF,SigT,Z,eGT,EcGM,EcRPA)
|
||||
|
||||
! Print one-electron energies and other stuff for G0T0
|
||||
|
||||
@ -8,7 +8,7 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: EGHF
|
||||
double precision,intent(in) :: EcGM
|
||||
double precision,intent(in) :: EcRPA
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
@ -48,14 +48,14 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
|
||||
end do
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO-LUMO gap = ',Gap*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO-LUMO gap = ',Gap*HaToeV,' eV'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy = ',EcRPA,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF total energy = ',ENuc + ERHF + EcRPA,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF correlation energy = ',EcGM,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF total energy = ',ENuc + ERHF + EcGM,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF correlation energy = ',EcRPA,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF total energy = ',ENuc + EGHF + EcRPA,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF correlation energy = ',EcGM,' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF total energy = ',ENuc + EGHF + EcGM,' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
||||
@ -33,7 +33,6 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
logical :: print_T = .false.
|
||||
logical :: dRPA
|
||||
integer :: ispin
|
||||
integer :: iblock
|
||||
integer :: nOOs,nOOt
|
||||
integer :: nVVs,nVVt
|
||||
double precision :: EcRPA(nspin)
|
||||
@ -74,18 +73,12 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
! Dimensions of the ppRPA linear reponse matrices
|
||||
|
||||
! nOOs = nO*(nO + 1)/2
|
||||
! nVVs = nV*(nV + 1)/2
|
||||
|
||||
nOOs = nO*nO
|
||||
nVVs = nV*nV
|
||||
nOOs = nO*(nO + 1)/2
|
||||
nVVs = nV*(nV + 1)/2
|
||||
|
||||
nOOt = nO*(nO - 1)/2
|
||||
nVVt = nV*(nV - 1)/2
|
||||
|
||||
! nOO = nO*nO
|
||||
! nVV = nV*nV
|
||||
|
||||
! Dimension of the supermatrix
|
||||
|
||||
n2h1p = (nOOs+nOOt)*nV
|
||||
@ -117,16 +110,14 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
! alpha-beta block
|
||||
|
||||
ispin = 1
|
||||
! iblock = 1
|
||||
iblock = 3
|
||||
|
||||
! Compute linear response
|
||||
|
||||
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
|
||||
|
||||
call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
|
||||
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
|
||||
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,eHF,ERI,Cpp)
|
||||
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
|
||||
|
||||
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
|
||||
|
||||
@ -135,23 +126,21 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
! Compute excitation densities
|
||||
|
||||
call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
call RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp,X1s,Y1s,X2s,Y2s)
|
||||
|
||||
! alpha-alpha block
|
||||
|
||||
ispin = 2
|
||||
! iblock = 2
|
||||
iblock = 4
|
||||
|
||||
! Compute linear response
|
||||
|
||||
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
|
||||
|
||||
call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
|
||||
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
|
||||
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
|
||||
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,eHF,ERI,Cpp)
|
||||
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
|
||||
|
||||
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
|
||||
|
||||
@ -160,10 +149,14 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
! Compute excitation densities
|
||||
|
||||
call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
call RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp,X1t,Y1t,X2t,Y2t)
|
||||
|
||||
else
|
||||
|
||||
allocate(rho1s(0,0,0),rho1t(0,0,0),rho2s(0,0,0),rho2t(0,0,0))
|
||||
|
||||
end if
|
||||
|
||||
! Memory allocation
|
||||
|
||||
@ -51,28 +51,23 @@ subroutine GGW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(b)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
num = (rho(a,i,m)*rho(b,j,m) - rho(b,i,m)*rho(a,j,m))/2
|
||||
dem = - Om(m) - eGW(b) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(a)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(a)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(a) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(b)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
dem = - Om(m) - eGW(a) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
dem = - Om(m) - eGW(b) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
|
||||
KB_dyn(ab,ij) = 0.5d0*KB_dyn(ab,ij)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
@ -53,36 +53,28 @@ subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om,
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = OmBSE - (eGW(a) + eGW(c) + Om(m))
|
||||
! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
|
||||
num = - rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
num = (rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))/2
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - (eGW(b) + eGW(d) + Om(m))
|
||||
! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
|
||||
num = - rho(b,c,m)*rho(a,d,m)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - (eGW(b) + eGW(c) + Om(m))
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - (eGW(a) + eGW(d) + Om(m))
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
|
||||
KC_dyn(ab,cd) = 0.5d0*KC_dyn(ab,cd)
|
||||
ZC_dyn(ab,cd) = 0.5d0*ZC_dyn(ab,cd)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
@ -100,7 +100,7 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
|
||||
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
|
||||
|
||||
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
|
||||
call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
|
||||
|
||||
! Singlet manifold
|
||||
|
||||
@ -130,7 +130,7 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
|
||||
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
|
||||
|
||||
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
|
||||
call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
|
||||
|
||||
end if
|
||||
|
||||
@ -187,7 +187,7 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
|
||||
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
|
||||
|
||||
call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
|
||||
call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
|
||||
|
||||
end if
|
||||
|
||||
|
||||
@ -90,7 +90,7 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
|
||||
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
|
||||
|
||||
call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
|
||||
call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
|
||||
if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
|
||||
|
||||
!-------------------
|
||||
! Singlet manifold
|
||||
@ -116,8 +116,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
|
||||
write(*,*)
|
||||
|
||||
call RGW_phBSE_static_kernel(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W)
|
||||
call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta)
|
||||
|
||||
call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta)
|
||||
if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta)
|
||||
|
||||
deallocate(W)
|
||||
|
||||
@ -137,55 +137,23 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
|
||||
allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
|
||||
|
||||
!print*, 'RGW_ppBSE_static_kernel_C:'
|
||||
!call wall_time(tt1)
|
||||
call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_C (sec)', tt2-tt1
|
||||
|
||||
!print*, 'RGW_ppBSE_static_kernel_D:'
|
||||
!call wall_time(tt1)
|
||||
call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_D (sec)', tt2-tt1
|
||||
|
||||
if(.not.TDA) then
|
||||
!print*, 'RGW_ppBSE_static_kernel_B:'
|
||||
!call wall_time(tt1)
|
||||
call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_B (sec)', tt2-tt1
|
||||
endif
|
||||
|
||||
!print*, 'ppLR_C:'
|
||||
!call wall_time(tt1)
|
||||
call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for ppLR_C (sec)', tt2-tt1
|
||||
|
||||
!print*, 'ppLR_D:'
|
||||
!call wall_time(tt1)
|
||||
call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for ppLR_D (sec)', tt2-tt1
|
||||
|
||||
if(.not.TDA) then
|
||||
!print*, 'ppLR_B:'
|
||||
!call wall_time(tt1)
|
||||
call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for ppLR_B (sec)', tt2-tt1
|
||||
endif
|
||||
|
||||
Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
|
||||
Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
|
||||
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
|
||||
|
||||
!print*, 'ppLR:'
|
||||
!call wall_time(tt1)
|
||||
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A,1X,F10.3)'), 'wall time for ppLR (sec)', tt2-tt1
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
!
|
||||
@ -197,18 +165,18 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
|
||||
|
||||
|
||||
! ---
|
||||
! Davidson
|
||||
! ---
|
||||
!! ---
|
||||
!! Davidson
|
||||
!! ---
|
||||
|
||||
!n_states = nOO + 5
|
||||
!n_states_diag = n_states + 4
|
||||
!allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag))
|
||||
!
|
||||
! supp_data_int = 1
|
||||
|
||||
!supp_data_int_size = 1
|
||||
!allocate(supp_data_int(supp_data_int_size))
|
||||
!supp_data_int(1) = nS
|
||||
!
|
||||
|
||||
!supp_data_dbl_size = nS + nOrb*nOrb*nS + 1
|
||||
!allocate(supp_data_dbl(supp_data_dbl_size))
|
||||
!! scalars
|
||||
@ -228,7 +196,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
! i_data = i_data + 1
|
||||
! supp_data_dbl(i_data) = OmRPA(m)
|
||||
!enddo
|
||||
!
|
||||
|
||||
!call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, &
|
||||
! 1.d0, & ! lambda
|
||||
! eGW(1), &
|
||||
@ -236,8 +204,8 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
! ERI(1,1,1,1), &
|
||||
! supp_data_int(1), supp_data_int_size, &
|
||||
! supp_data_dbl(1), supp_data_dbl_size, &
|
||||
! Om(1), R(1,1), n_states, n_states_diag, "GW")
|
||||
!
|
||||
! Om(1), R(1,1), n_states, n_states_diag, "GW", 1)
|
||||
|
||||
!deallocate(Om, R)
|
||||
!deallocate(supp_data_dbl, supp_data_int)
|
||||
!stop
|
||||
@ -308,10 +276,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
!print*, 'LAPACK:'
|
||||
!print*, Om2
|
||||
!print*, Om1
|
||||
|
||||
! ---
|
||||
! Davidson
|
||||
! ---
|
||||
@ -320,7 +284,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
!n_states_diag = n_states + 4
|
||||
!allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag))
|
||||
|
||||
!supp_data_int = 1
|
||||
!supp_data_int_size = 1
|
||||
!allocate(supp_data_int(supp_data_int_size))
|
||||
!supp_data_int(1) = nS
|
||||
|
||||
@ -351,7 +315,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
! ERI(1,1,1,1), &
|
||||
! supp_data_int(1), supp_data_int_size, &
|
||||
! supp_data_dbl(1), supp_data_dbl_size, &
|
||||
! Om(1), R(1,1), n_states, n_states_diag, "GW")
|
||||
! Om(1), R(1,1), n_states, n_states_diag, "GW", 1)
|
||||
|
||||
!deallocate(Om, R)
|
||||
!deallocate(supp_data_dbl, supp_data_int)
|
||||
@ -371,15 +335,9 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
|
||||
call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
|
||||
Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
|
||||
|
||||
|
||||
!----------------!
|
||||
! Upfolded ppBSE !
|
||||
!----------------!
|
||||
|
||||
! call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
|
||||
|
||||
deallocate(KB_sta,KC_sta,KD_sta)
|
||||
deallocate(Om1,X1,Y1,Om2,X2,Y2)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -53,30 +53,23 @@ subroutine RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lamb
|
||||
ij = ij + 1
|
||||
|
||||
do m=1,nS
|
||||
num = (rho(a,i,m)*rho(b,j,m) + rho(b,i,m)*rho(a,j,m))/2
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(b)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(b) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(a)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(a) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(a)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(a) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(b)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(b) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
|
||||
KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -99,31 +92,22 @@ subroutine RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lamb
|
||||
ij = ij + 1
|
||||
|
||||
do m=1,nS
|
||||
num = (rho(a,i,m)*rho(b,j,m) - rho(b,i,m)*rho(a,j,m))/2
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(b)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(b) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(a)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(a)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
dem = - Om(m) - eGW(a) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(b)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
dem = - Om(m) - eGW(a) + eGW(j)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
dem = - Om(m) - eGW(b) + eGW(i)
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
|
||||
KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
@ -55,35 +55,28 @@ subroutine RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,e
|
||||
cd = cd + 1
|
||||
|
||||
do m=1,nS
|
||||
num = (rho(a,c,m)*rho(b,d,m) + rho(b,c,m)*rho(a,d,m))/2
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(b)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(a)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(a)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(b)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
|
||||
KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -106,36 +99,26 @@ subroutine RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,e
|
||||
cd = cd + 1
|
||||
|
||||
do m=1,nS
|
||||
num = (rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))/2
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(b)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(a)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(a)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(b)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
dem = OmBSE - Om(m) - eGW(a) - eGW(d)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
dem = OmBSE - Om(m) - eGW(b) - eGW(c)
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
|
||||
KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)
|
||||
ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
@ -57,7 +57,6 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
|
||||
do m=1,nS
|
||||
num = (rho(i,k,m)*rho(j,l,m) + rho(j,k,m)*rho(i,l,m))/2
|
||||
|
||||
! dem = - Om(m)
|
||||
dem = - OmBSE - Om(m) + eGW(j) + eGW(l)
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
@ -101,6 +100,7 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
|
||||
|
||||
do m=1,nS
|
||||
num = (rho(i,k,m)*rho(j,l,m) - rho(j,k,m)*rho(i,l,m))/2
|
||||
|
||||
dem = - OmBSE - Om(m) + eGW(j) + eGW(l)
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
@ -119,9 +119,6 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
|
||||
|
||||
end do
|
||||
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
@ -1,263 +0,0 @@
|
||||
subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
|
||||
|
||||
! Upfolded ppBSE@GW (TDA only!)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: rho(nOrb,nOrb,nS)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: eGW(nOrb)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: s
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
integer :: m,ij,kl,ijm
|
||||
integer,parameter :: maxH = 100
|
||||
double precision :: tmp,tmp1,tmp2,tmp3,tmp4
|
||||
|
||||
integer :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,n2h2p,nH
|
||||
double precision,external :: Kronecker_delta
|
||||
|
||||
integer,allocatable :: order(:)
|
||||
double precision,allocatable :: H(:,:)
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: OmBSE(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)'* Upfolded ppBSE@GW Calculation *'
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)
|
||||
|
||||
! TDA for W
|
||||
|
||||
write(*,*) 'Tamm-Dancoff approximation by default!'
|
||||
write(*,*)
|
||||
|
||||
! Dimension of the supermatrix
|
||||
|
||||
n1h = nO
|
||||
n1p = nV
|
||||
|
||||
if(ispin == 1) then
|
||||
|
||||
n2h = nO*(nO+1)/2
|
||||
n2p = nV*(nV+1)/2
|
||||
|
||||
end if
|
||||
|
||||
if(ispin == 2) then
|
||||
|
||||
n2h = nO*(nO-1)/2
|
||||
n2p = nV*(nV-1)/2
|
||||
|
||||
end if
|
||||
|
||||
n1h1p = n1h*n1p
|
||||
|
||||
n3h1p = n2h*n1h1p
|
||||
n3p1h = n2p*n1h1p
|
||||
|
||||
nH = n2h + 4*n3h1p
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH))
|
||||
|
||||
! Initialization
|
||||
|
||||
H(:,:) = 0d0
|
||||
|
||||
!----------------------------------------!
|
||||
! Compute BSE supermatrix !
|
||||
!----------------------------------------!
|
||||
! !
|
||||
! | D -M1 -M1 -M2 -M2 | !
|
||||
! | | !
|
||||
! | +M1 E1 0 0 0 | !
|
||||
! | | !
|
||||
! H = | +M1 0 E2 0 0 | !
|
||||
! | | !
|
||||
! | +M2 0 0 E3 0 | !
|
||||
! | | !
|
||||
! | +M2 0 0 0 E4 | !
|
||||
! !
|
||||
!----------------------------------------!
|
||||
|
||||
!----------------------!
|
||||
! Block D for singlets !
|
||||
!----------------------!
|
||||
|
||||
if(ispin == 1) then
|
||||
|
||||
ij = 0
|
||||
do i=nC+1,nO
|
||||
do j=i,nO
|
||||
ij = ij + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k,nO
|
||||
kl = kl + 1
|
||||
|
||||
H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
|
||||
+ (ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
!----------------------!
|
||||
! Block D for triplets !
|
||||
!----------------------!
|
||||
|
||||
if(ispin == 2) then
|
||||
|
||||
ij = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
|
||||
+ (ERI(i,j,k,l) - ERI(i,j,l,k))
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
!----------------!
|
||||
! Blocks M1 & M2 !
|
||||
!----------------!
|
||||
|
||||
ijm = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
do m=1,nS
|
||||
ijm = ijm + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
tmp1 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,k,m)
|
||||
tmp2 = sqrt(2d0)*Kronecker_delta(j,k)*rho(i,l,m)
|
||||
tmp3 = sqrt(2d0)*Kronecker_delta(i,l)*rho(j,k,m)
|
||||
tmp4 = sqrt(2d0)*Kronecker_delta(i,k)*rho(j,l,m)
|
||||
|
||||
H(n2h+0*n3h1p+ijm,kl ) = +tmp1
|
||||
H(kl ,n2h+0*n3h1p+ijm) = +tmp4
|
||||
|
||||
H(n2h+1*n3h1p+ijm,kl ) = +tmp1
|
||||
H(kl ,n2h+1*n3h1p+ijm) = +tmp4
|
||||
|
||||
H(n2h+2*n3h1p+ijm,kl ) = +tmp2
|
||||
H(kl ,n2h+2*n3h1p+ijm) = +tmp3
|
||||
|
||||
H(n2h+3*n3h1p+ijm,kl ) = +tmp2
|
||||
H(kl ,n2h+3*n3h1p+ijm) = +tmp4
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!------------!
|
||||
! Block 3h1p !
|
||||
!------------!
|
||||
|
||||
ijm = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
do m=1,nS
|
||||
ijm = ijm + 1
|
||||
|
||||
tmp = - eGW(i) - eGW(j) + Om(m)
|
||||
|
||||
H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp
|
||||
H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp
|
||||
H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp
|
||||
H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------------------!
|
||||
! Diagonalize supermatrix !
|
||||
!-------------------------!
|
||||
|
||||
call diagonalize_general_matrix(nH,H,OmBSE,X)
|
||||
|
||||
do s=1,nH
|
||||
order(s) = s
|
||||
end do
|
||||
|
||||
call quick_sort(OmBSE,order,nH)
|
||||
call set_order(X,order,nH,nH)
|
||||
|
||||
!-----------------!
|
||||
! Compute weights !
|
||||
!-----------------!
|
||||
|
||||
Z(:) = 0d0
|
||||
do s=1,nH
|
||||
do ij=1,n2h
|
||||
Z(s) = Z(s) + X(ij,s)**2
|
||||
end do
|
||||
end do
|
||||
|
||||
!--------------!
|
||||
! Dump results !
|
||||
!--------------!
|
||||
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)' Upfolded ppBSE excitation energies (eV) '
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
|
||||
'|','#','|','OmBSE (eV)','|','Z','|'
|
||||
write(*,*)'-------------------------------------------'
|
||||
|
||||
do s=1,min(nH,maxH)
|
||||
if(Z(s) > 1d-7) &
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',s,'|',-OmBSE(s)*HaToeV,'|',Z(s),'|'
|
||||
end do
|
||||
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
@ -194,10 +194,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
@ -210,10 +210,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
|
||||
write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
||||
@ -113,13 +113,13 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,
|
||||
allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc))
|
||||
|
||||
call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA)
|
||||
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB)
|
||||
|
||||
! Spin-conserved manifold
|
||||
|
||||
@ -145,23 +145,23 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,
|
||||
if(doXBS) then
|
||||
|
||||
call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA)
|
||||
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
|
||||
|
||||
end if
|
||||
|
||||
call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
Aph(:,:) = Aph(:,:) + KA(:,:)
|
||||
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
|
||||
|
||||
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
|
||||
call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
|
||||
|
||||
call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
|
||||
@ -209,23 +209,23 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,
|
||||
if(doXBS) then
|
||||
|
||||
call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA)
|
||||
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
|
||||
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
|
||||
|
||||
end if
|
||||
|
||||
call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
Aph(:,:) = Aph(:,:) + KA(:,:)
|
||||
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
|
||||
|
||||
call phULR(TDA_W,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
|
||||
call phULR(TDA,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
|
||||
|
||||
call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
|
||||
|
||||
@ -121,7 +121,7 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
|
||||
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB)
|
||||
|
||||
Aph(:,:) = Aph(:,:) + KA(:,:)
|
||||
! if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
|
||||
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
|
||||
|
||||
call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
|
||||
|
||||
@ -155,6 +155,7 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
|
||||
allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf))
|
||||
|
||||
! Compute spin-conserved BSE excitation energies
|
||||
|
||||
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
|
||||
@ -32,6 +32,10 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
|
||||
double precision,intent(out) :: KA(nSt,nSt)
|
||||
|
||||
! Initialization
|
||||
|
||||
KA(:,:) = 0d0
|
||||
|
||||
!--------------------------------------------------!
|
||||
! Build BSE matrix for spin-conserving transitions !
|
||||
!--------------------------------------------------!
|
||||
@ -55,7 +59,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi
|
||||
KA(ia,jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -79,7 +83,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi
|
||||
KA(nSa+ia,nSa+jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -111,7 +115,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi
|
||||
KA(ia,jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -135,7 +139,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi
|
||||
KA(nSa+ia,nSa+jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -32,6 +32,10 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
|
||||
double precision,intent(out) :: KB(nSt,nSt)
|
||||
|
||||
! Initialization
|
||||
|
||||
KB(:,:) = 0d0
|
||||
|
||||
!--------------------------------------------------!
|
||||
! Build BSE matrix for spin-conserving transitions !
|
||||
!--------------------------------------------------!
|
||||
@ -55,7 +59,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KB(ia,jb) = KB(ia,jb) + 2d0*lambda*chi
|
||||
KB(ia,jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -80,7 +84,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KB(nSa+ia,nSa+jb) = KB(nSa+ia,nSa+jb) + 2d0*lambda*chi
|
||||
KB(nSa+ia,nSa+jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -113,7 +117,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KB(ia,nSa+jb) = KB(ia,nSa+jb) + 2d0*lambda*chi
|
||||
KB(ia,nSa+jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -137,7 +141,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
|
||||
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Om(kc)/eps
|
||||
end do
|
||||
|
||||
KB(nSa+ia,jb) = KB(nSa+ia,jb) + 2d0*lambda*chi
|
||||
KB(nSa+ia,jb) = 2d0*lambda**2*chi
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -85,7 +85,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
|
||||
|
||||
! SRG regularization
|
||||
|
||||
flow = 500d0
|
||||
flow = 100d0
|
||||
|
||||
if(doSRG) then
|
||||
|
||||
|
||||
@ -111,6 +111,10 @@ subroutine ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
deallocate(Aph,Bph,XpY,XmY)
|
||||
|
||||
else
|
||||
|
||||
allocate(rho(0,0,0))
|
||||
|
||||
end if
|
||||
|
||||
!-------------------------!
|
||||
|
||||
@ -111,6 +111,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
deallocate(Aph,Bph,XpY,XmY)
|
||||
|
||||
else
|
||||
|
||||
allocate(rho(0,0,0))
|
||||
|
||||
end if
|
||||
|
||||
! Initialization
|
||||
|
||||
@ -31,7 +31,6 @@ subroutine print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
|
||||
double precision :: Sz
|
||||
double precision :: Sx2,Sy2,Sz2
|
||||
integer :: mu,nu
|
||||
double precision,allocatable :: qa(:),qb(:)
|
||||
|
||||
logical :: dump_orb = .false.
|
||||
|
||||
@ -119,18 +118,4 @@ subroutine print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
|
||||
call vecout(nBas,eHF(:,2))
|
||||
write(*,*)
|
||||
|
||||
allocate(qa(nBas),qb(nBas))
|
||||
|
||||
qa(:) = 0d0
|
||||
qb(:) = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
qa(mu) = qa(mu) + P(mu,nu,1)*S(nu,mu)
|
||||
qb(mu) = qb(mu) + P(mu,nu,2)*S(nu,mu)
|
||||
end do
|
||||
end do
|
||||
|
||||
call vecout(nBas,qa)
|
||||
call vecout(nBas,qb)
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -39,6 +39,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
|
||||
delta_dRPA = 0d0
|
||||
if(dRPA) delta_dRPA = 1d0
|
||||
|
||||
! Initialization
|
||||
|
||||
Aph(:,:) = 0d0
|
||||
|
||||
!----------------------------------------------
|
||||
! Build A matrix for spin-conserved transitions
|
||||
!----------------------------------------------
|
||||
@ -127,8 +131,6 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E
|
||||
|
||||
if(ispin == 2) then
|
||||
|
||||
Aph(:,:) = 0d0
|
||||
|
||||
! abab block
|
||||
|
||||
ia = 0
|
||||
|
||||
@ -38,6 +38,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
|
||||
delta_dRPA = 0d0
|
||||
if(dRPA) delta_dRPA = 1d0
|
||||
|
||||
! Initialization
|
||||
|
||||
Bph(:,:) = 0d0
|
||||
|
||||
!-----------------------------------------------
|
||||
! Build B matrix for spin-conserving transitions
|
||||
!-----------------------------------------------
|
||||
@ -124,8 +128,6 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a
|
||||
|
||||
if(ispin == 2) then
|
||||
|
||||
Bph(:,:) = 0d0
|
||||
|
||||
! abba block
|
||||
|
||||
ia = 0
|
||||
|
||||
@ -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))
|
||||
|
||||
! Hermitian case for TDA
|
||||
|
||||
if(TDA) then
|
||||
|
||||
X1(:,:) = +Cpp(:,:)
|
||||
@ -89,8 +88,8 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
|
||||
! imp_bio = .True. ! impose bi-orthogonality
|
||||
! 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
|
||||
@ -100,7 +99,7 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
! Y2(j,i) = Z(nVV+j,i)
|
||||
! enddo
|
||||
! enddo
|
||||
!
|
||||
|
||||
! do i = 1, nVV
|
||||
! Om1(i) = Om(nOO+i)
|
||||
! do j = 1, nVV
|
||||
|
||||
@ -89,8 +89,8 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
|
||||
! imp_bio = .True. ! impose bi-orthogonality
|
||||
! 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
|
||||
@ -100,7 +100,7 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
! Y2(j,i) = Z(nVV+j,i)
|
||||
! enddo
|
||||
! enddo
|
||||
!
|
||||
|
||||
! do i = 1, nVV
|
||||
! Om1(i) = Om(nOO+i)
|
||||
! do j = 1, nVV
|
||||
|
||||
@ -33,7 +33,7 @@ subroutine ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF,
|
||||
|
||||
call wall_time(t1)
|
||||
|
||||
if((nOO+nVV) .le. 20000) then
|
||||
if((nOO+nVV) .le. 30000) then
|
||||
|
||||
call ppLR_GW_HR_calc_oneshot(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, n_states_diag, &
|
||||
ERI(1,1,1,1), eta, rho(1,1,1), Om(1), U(1,1), W(1,1))
|
||||
|
||||
@ -4,7 +4,7 @@
|
||||
subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, &
|
||||
supp_data_int, supp_data_int_size, &
|
||||
supp_data_dbl, supp_data_dbl_size, &
|
||||
Om, R, n_states, n_states_diag, kernel)
|
||||
Om, R, n_states, n_states_diag, kernel, mode_dav)
|
||||
|
||||
!
|
||||
! Extract the low n_states
|
||||
@ -18,6 +18,54 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
! (-B.T -D)
|
||||
!
|
||||
|
||||
implicit none
|
||||
|
||||
logical, intent(in) :: TDA
|
||||
integer, intent(in) :: ispin
|
||||
integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV
|
||||
integer, intent(in) :: n_states ! nb of physical states
|
||||
integer, intent(in) :: n_states_diag ! nb of states used to get n_states
|
||||
integer, intent(in) :: supp_data_int_size
|
||||
integer, intent(in) :: supp_data_dbl_size
|
||||
integer, intent(in) :: mode_dav
|
||||
character(len=*), intent(in) :: kernel
|
||||
double precision, intent(in) :: lambda, eF
|
||||
double precision, intent(in) :: e(nOrb)
|
||||
double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
integer, intent(in) :: supp_data_int(supp_data_int_size)
|
||||
double precision, intent(in) :: supp_data_dbl(supp_data_dbl_size)
|
||||
double precision, intent(out) :: Om(n_states)
|
||||
double precision, intent(out) :: R(nOO+nVV,n_states_diag)
|
||||
|
||||
if(mode_dav .eq. 1) then
|
||||
|
||||
call ppLR_davidson_1(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e(1), eF, ERI(1,1,1,1), &
|
||||
supp_data_int(1), supp_data_int_size, supp_data_dbl(1), supp_data_dbl_size, &
|
||||
Om(1), R(1,1), n_states, n_states_diag, kernel)
|
||||
|
||||
elseif(mode_dav .eq. 2) then
|
||||
|
||||
call ppLR_davidson_2(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e(1), eF, ERI(1,1,1,1), &
|
||||
supp_data_int(1), supp_data_int_size, supp_data_dbl(1), supp_data_dbl_size, &
|
||||
Om(1), R(1,1), n_states, n_states_diag, kernel)
|
||||
|
||||
else
|
||||
|
||||
print*, " unknown Davidson's variant"
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine ppLR_davidson_1(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, &
|
||||
supp_data_int, supp_data_int_size, &
|
||||
supp_data_dbl, supp_data_dbl_size, &
|
||||
Om, R, n_states, n_states_diag, kernel)
|
||||
|
||||
use omp_lib
|
||||
|
||||
implicit none
|
||||
@ -43,16 +91,15 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
integer :: shift1, shift2
|
||||
integer :: i, j, k, l, ab
|
||||
integer :: p, q, mm, i_data, nS
|
||||
integer :: i_omax(n_states)
|
||||
logical :: converged
|
||||
character*(16384) :: write_buffer
|
||||
character(len=6+41*n_states) :: write_buffer
|
||||
double precision :: r1, r2, dtwo_pi
|
||||
double precision :: lambda_tmp
|
||||
double precision :: to_print(2,n_states)
|
||||
double precision :: mem
|
||||
double precision :: eta
|
||||
double precision :: t1, t2, tt1, tt2
|
||||
character(len=len(kernel)) :: kernel_name
|
||||
integer, allocatable :: i_omax(:)
|
||||
double precision, allocatable :: H_diag(:)
|
||||
double precision, allocatable :: W(:,:)
|
||||
double precision, allocatable :: U(:,:)
|
||||
@ -61,6 +108,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
double precision, allocatable :: overlap(:)
|
||||
double precision, allocatable :: S_check(:,:)
|
||||
double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:)
|
||||
double precision, allocatable :: to_print(:,:)
|
||||
|
||||
double precision, external :: u_dot_u
|
||||
|
||||
@ -92,7 +140,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
write(*,'(A40, A12)') 'Kernel: ', kernel_name
|
||||
|
||||
|
||||
|
||||
allocate(i_omax(n_states))
|
||||
allocate(to_print(2,n_states))
|
||||
allocate(H_diag(N))
|
||||
allocate(U(N,M))
|
||||
allocate(W(N,M))
|
||||
@ -101,7 +150,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
allocate(residual_norm(n_states_diag))
|
||||
|
||||
mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) &
|
||||
+ 8.d0 * dble(supp_data_dbl_size) + 4.d0 * dble(supp_data_int_size)
|
||||
+ 8.d0 * dble(2*supp_data_dbl_size) + 4.d0 * dble(2*supp_data_int_size)
|
||||
|
||||
write(*,'(A40, F12.4)') 'I/O mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0)
|
||||
|
||||
@ -118,6 +167,10 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, &
|
||||
ERI(1,1,1,1), H_diag(1))
|
||||
|
||||
! to avoid compiler warnings
|
||||
allocate(rho_tmp(0,0,0))
|
||||
allocate(Om_tmp(0))
|
||||
|
||||
elseif(kernel_name .eq. "gw") then
|
||||
|
||||
nS = supp_data_int(1)
|
||||
@ -214,6 +267,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
write(6,'(A)') write_buffer(1:6+41*n_states)
|
||||
|
||||
|
||||
W = 0.d0
|
||||
converged = .False.
|
||||
itertot = 0
|
||||
|
||||
@ -325,7 +379,6 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, W(1,shift2+1), size(W, 1))
|
||||
|
||||
! check if W1 = U1 h_val
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i, k) &
|
||||
!$OMP SHARED(n_states, n_states_diag, N, shift2, U, h_val, W, H_diag, residual_norm, to_print)
|
||||
@ -342,13 +395,11 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!print*, " to_print", to_print
|
||||
|
||||
if((itertot > 1) .and. (iter == 1)) then
|
||||
continue
|
||||
else
|
||||
write(*,'(1X, I3, 1X, 100(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states)
|
||||
!write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, F16.10, 1X, F16.10))') iter-1, to_print(1:2,1:n_states)
|
||||
write(*,'(1X, I3, 1X, 10000(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states)
|
||||
endif
|
||||
|
||||
!call wall_time(tt2)
|
||||
@ -361,7 +412,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
endif
|
||||
|
||||
do k = 1, n_states
|
||||
if(residual_norm(k) > 1.d8) then
|
||||
if(residual_norm(k) > 1.d10) then
|
||||
print *, 'Davidson failed'
|
||||
stop -1
|
||||
endif
|
||||
@ -377,7 +428,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
|
||||
call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, &
|
||||
W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, R, size(R, 1))
|
||||
0.d0, R(1,1), size(R, 1))
|
||||
|
||||
do k = 1, n_states_diag
|
||||
do i = 1, N
|
||||
W(i,k) = R(i,k)
|
||||
@ -427,6 +479,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
|
||||
print*, k, Om(k)
|
||||
enddo
|
||||
|
||||
deallocate(i_omax)
|
||||
deallocate(to_print)
|
||||
deallocate(H_diag)
|
||||
deallocate(U)
|
||||
deallocate(W)
|
||||
@ -449,3 +503,404 @@ end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine ppLR_davidson_2(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, &
|
||||
supp_data_int, supp_data_int_size, &
|
||||
supp_data_dbl, supp_data_dbl_size, &
|
||||
Om, R, n_states, n_states_diag, kernel)
|
||||
|
||||
use omp_lib
|
||||
|
||||
implicit none
|
||||
|
||||
logical, intent(in) :: TDA
|
||||
integer, intent(in) :: ispin
|
||||
integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV
|
||||
integer, intent(in) :: n_states ! nb of physical states
|
||||
integer, intent(in) :: n_states_diag ! nb of states used to get n_states
|
||||
integer, intent(in) :: supp_data_int_size
|
||||
integer, intent(in) :: supp_data_dbl_size
|
||||
character(len=*), intent(in) :: kernel
|
||||
double precision, intent(in) :: lambda, eF
|
||||
double precision, intent(in) :: e(nOrb)
|
||||
double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
integer, intent(in) :: supp_data_int(supp_data_int_size)
|
||||
double precision, intent(in) :: supp_data_dbl(supp_data_dbl_size)
|
||||
double precision, intent(out) :: Om(n_states)
|
||||
double precision, intent(out) :: R(nOO+nVV,n_states_diag)
|
||||
|
||||
integer :: N, M, num_threads, n_states_delta
|
||||
integer :: it_start, it_delta, it_size
|
||||
integer :: iter, itermax, itertot
|
||||
integer :: i, j, k, l, ab
|
||||
integer :: p, q, mm, i_data, nS
|
||||
logical :: converged
|
||||
double precision :: r1, r2, dtwo_pi
|
||||
double precision :: mem
|
||||
double precision :: eta
|
||||
double precision :: t1, t2, tt1, tt2
|
||||
character(len=len(kernel)) :: kernel_name
|
||||
integer, allocatable :: i_omax(:)
|
||||
character(len=:), allocatable :: write_buffer
|
||||
double precision, allocatable :: to_print(:,:)
|
||||
double precision, allocatable :: H_diag(:)
|
||||
double precision, allocatable :: W0(:,:), W1(:,:)
|
||||
double precision, allocatable :: U0(:,:), U1(:,:)
|
||||
double precision, allocatable :: h(:,:), h_vec(:,:), h_val(:)
|
||||
double precision, allocatable :: residual_norm(:)
|
||||
double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:)
|
||||
|
||||
double precision, external :: u_dot_u
|
||||
|
||||
call wall_time(t1)
|
||||
|
||||
dtwo_pi = 6.283185307179586d0
|
||||
|
||||
N = nOO + nVV
|
||||
|
||||
n_states_delta = min(max(25, n_states_diag/2), n_states_diag)
|
||||
itermax = 8
|
||||
M = n_states_diag + itermax * n_states_delta
|
||||
|
||||
call lower_case(trim(kernel), kernel_name)
|
||||
|
||||
if(M .ge. N) then
|
||||
print*, 'N = ', N
|
||||
print*, 'M = ', M
|
||||
print*, ' use Lapack or decrease n_states and/or itermax '
|
||||
stop
|
||||
endif
|
||||
|
||||
write(6,'(A)') ''
|
||||
write(6,'(A)') 'Davidson Diagonalization'
|
||||
write(6,'(A)') '------------------------'
|
||||
write(6,'(A)') ''
|
||||
|
||||
write(*,'(A40, I12)') 'Number of states = ', n_states
|
||||
write(*,'(A40, I12)') 'Number of states in diag = ', n_states_diag
|
||||
write(*,'(A40, I12)') 'Number of states to add = ', n_states_delta
|
||||
write(*,'(A40, I12)') 'Number of basis functions = ', N
|
||||
write(*,'(A40, A12)') 'Kernel: ', kernel_name
|
||||
|
||||
|
||||
|
||||
allocate(character(len=50*n_states) :: write_buffer)
|
||||
allocate(i_omax(n_states))
|
||||
allocate(to_print(2,n_states))
|
||||
allocate(H_diag(N))
|
||||
allocate(U0(N,M), U1(N,n_states_diag))
|
||||
allocate(W0(N,M), W1(N,n_states_diag))
|
||||
allocate(h(M,M), h_vec(M,M), h_val(M))
|
||||
allocate(residual_norm(n_states_diag))
|
||||
|
||||
mem = 8.d0 * dble(nOrb) + 8.d0 * dble(nOrb)**4 + 8.d0 * dble(N*n_states) &
|
||||
+ 8.d0 * dble(2*supp_data_dbl_size) + 4.d0 * dble(2*supp_data_int_size)
|
||||
|
||||
write(*,'(A40, F12.4)') 'I/O mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0)
|
||||
|
||||
mem = 8.d0 * dble(N) &
|
||||
+ 8.d0 * dble(N*M) &
|
||||
+ 8.d0 * dble(N*M) &
|
||||
+ 8.d0 * dble(N*n_states_diag) &
|
||||
+ 8.d0 * dble(N*n_states_diag) &
|
||||
+ 8.d0 * dble(M*M) &
|
||||
+ 8.d0 * dble(M*M) &
|
||||
+ 8.d0 * dble(M) &
|
||||
+ 8.d0 * dble(n_states_diag) &
|
||||
+ 1.d0 * dble(50*n_states)
|
||||
|
||||
write(*,'(A40, F12.4)') 'tmp mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0)
|
||||
|
||||
num_threads = omp_get_max_threads()
|
||||
write(*,'(A40, I12)') 'Number of threads = ', num_threads
|
||||
|
||||
|
||||
if(kernel_name .eq. "rpa") then
|
||||
|
||||
allocate(rho_tmp(0,0,0))
|
||||
allocate(Om_tmp(0))
|
||||
|
||||
call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, &
|
||||
ERI(1,1,1,1), H_diag(1))
|
||||
|
||||
elseif(kernel_name .eq. "gw") then
|
||||
|
||||
nS = supp_data_int(1)
|
||||
|
||||
allocate(rho_tmp(nS,nOrb,nOrb))
|
||||
allocate(Om_tmp(nS))
|
||||
|
||||
eta = supp_data_dbl(1)
|
||||
i_data = 1
|
||||
do q = 1, nOrb
|
||||
do p = 1, nOrb
|
||||
do mm = 1, nS
|
||||
i_data = i_data + 1
|
||||
rho_tmp(mm,p,q) = supp_data_dbl(i_data)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do mm = 1, nS
|
||||
i_data = i_data + 1
|
||||
Om_tmp(mm) = supp_data_dbl(i_data)
|
||||
enddo
|
||||
|
||||
call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, &
|
||||
ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), H_diag(1))
|
||||
|
||||
!! TODO
|
||||
!elseif(kernel_name .eq. "gf2") then
|
||||
|
||||
else
|
||||
|
||||
print*, ' kernel not supported', kernel
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
U0 = 0.d0
|
||||
W0 = 0.d0
|
||||
U1 = 0.d0
|
||||
W1 = 0.d0
|
||||
|
||||
! TODO: improve guess
|
||||
! initialize guess
|
||||
R = 0.d0
|
||||
do k = 1, n_states
|
||||
R(k,k) = 1.d0
|
||||
enddo
|
||||
do k = n_states+1, n_states_diag
|
||||
do i = 1, N
|
||||
call random_number(r1)
|
||||
call random_number(r2)
|
||||
r1 = dsqrt(-2.d0*dlog(r1))
|
||||
r2 = dtwo_pi*r2
|
||||
R(i,k) = r1*dcos(r2)
|
||||
enddo
|
||||
R(k,k) = R(k,k) + 10.d0
|
||||
call normalize(R(1,k), N)
|
||||
enddo
|
||||
|
||||
do k = 1, n_states_diag
|
||||
U0(:,k) = R(:,k)
|
||||
enddo
|
||||
|
||||
|
||||
write(6,'(A)') ''
|
||||
write_buffer = '====='
|
||||
do i = 1, n_states
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*n_states)
|
||||
write_buffer = 'Iter'
|
||||
do i = 1, n_states
|
||||
write_buffer = trim(write_buffer)//' Energy Residual '
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*n_states)
|
||||
write_buffer = '====='
|
||||
do i = 1, n_states
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') write_buffer(1:6+41*n_states)
|
||||
|
||||
|
||||
converged = .False.
|
||||
itertot = 0
|
||||
|
||||
do while (.not.converged)
|
||||
|
||||
itertot = itertot + 1
|
||||
if(itertot == itermax) then
|
||||
print*, 'exit before convergence !'
|
||||
print*, 'itertot == itermax', itertot
|
||||
exit
|
||||
endif
|
||||
|
||||
do iter = 1, itermax-1
|
||||
|
||||
if(iter .eq. 1) then
|
||||
it_start = 0
|
||||
it_delta = n_states_diag
|
||||
else
|
||||
it_start = n_states_diag + n_states_delta * (iter - 2)
|
||||
it_delta = n_states_delta
|
||||
endif
|
||||
|
||||
it_size = it_start + it_delta
|
||||
|
||||
if((iter > 1) .or. (itertot == 1)) then
|
||||
|
||||
!call wall_time(tt1)
|
||||
|
||||
call ortho_qr(U0(1,1), size(U0, 1), N, it_size)
|
||||
|
||||
if(kernel_name .eq. "rpa") then
|
||||
|
||||
call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, it_delta, &
|
||||
ERI(1,1,1,1), &
|
||||
U0(1,it_start+1), W0(1,it_start+1))
|
||||
|
||||
elseif(kernel_name .eq. "gw") then
|
||||
|
||||
call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, it_delta, &
|
||||
ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), &
|
||||
U0(1,it_start+1), W0(1,it_start+1))
|
||||
|
||||
!! TODO
|
||||
!elseif(kernel_name .eq. "gf2") then
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
! computed below
|
||||
continue
|
||||
endif
|
||||
|
||||
! h = U0.T H U0
|
||||
call dgemm('T', 'N', it_size, it_size, N, 1.d0, &
|
||||
U0(1,1), size(U0, 1), W0(1,1), size(W0, 1), &
|
||||
0.d0, h(1,1), size(h, 1))
|
||||
|
||||
! h h_vec = h_val h_vec
|
||||
call diag_nonsym_right(it_size, h(1,1), size(h, 1), h_vec(1,1), size(h_vec, 1), h_val(1), size(h_val, 1))
|
||||
|
||||
! U1(:,1:it_delta) = U0 h_vec(:,1:it_delta)
|
||||
call dgemm('N', 'N', N, it_delta, it_size, 1.d0, &
|
||||
U0(1,1), size(U0, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, U1(1,1), size(U1, 1))
|
||||
|
||||
do k = 1, it_delta
|
||||
call normalize(U1(1,k), N)
|
||||
enddo
|
||||
|
||||
! W1(:,1:it_delta) = W0 h_vec(:,1:it_delta)
|
||||
call dgemm('N', 'N', N, it_delta, it_size, 1.d0, &
|
||||
W0(1,1), size(W0, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, W1(1,1), size(W1, 1))
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i, k) &
|
||||
!$OMP SHARED(n_states, it_size, it_delta, N, U0, U1, &
|
||||
!$OMP h_val, W1, H_diag, residual_norm, to_print)
|
||||
!$OMP DO
|
||||
do k = 1, it_delta
|
||||
do i = 1, N
|
||||
U1(i,k) = (h_val(k) * U1(i,k) - W1(i,k)) / max(dabs(H_diag(i) - h_val(k)), 1.d-2)
|
||||
U0(i,it_size+k) = U1(i,k)
|
||||
enddo
|
||||
if(k <= n_states) then
|
||||
residual_norm(k) = u_dot_u(U1(1,k), N)
|
||||
to_print(1,k) = h_val(k)
|
||||
to_print(2,k) = residual_norm(k)
|
||||
endif
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
if((itertot > 1) .and. (iter == 1)) then
|
||||
continue
|
||||
else
|
||||
write(*,'(1X, I3, 1X, 10000(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states)
|
||||
endif
|
||||
|
||||
!call wall_time(tt2)
|
||||
!write(*,'(A50, F12.4)') 'wall time for one Davidson iteration (sec): ', tt2-tt1
|
||||
!stop
|
||||
|
||||
if(iter > 1) then
|
||||
converged = dabs(maxval(residual_norm(1:n_states))) < 1d-15
|
||||
endif
|
||||
|
||||
do k = 1, n_states
|
||||
if(residual_norm(k) > 1.d10) then
|
||||
print *, 'Davidson failed'
|
||||
stop -1
|
||||
endif
|
||||
enddo
|
||||
|
||||
if(converged) exit
|
||||
|
||||
enddo ! loop over iter
|
||||
|
||||
|
||||
! Re-contract U0 and update W0
|
||||
! --------------------------------
|
||||
|
||||
call dgemm('N', 'N', N, n_states_diag, it_size, 1.d0, &
|
||||
W0(1,1), size(W0, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, R(1,1), size(R, 1))
|
||||
|
||||
do k = 1, n_states_diag
|
||||
do i = 1, N
|
||||
W0(i,k) = R(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call dgemm('N', 'N', N, n_states_diag, it_size, 1.d0, &
|
||||
U0(1,1), size(U0, 1), h_vec(1,1), size(h_vec, 1), &
|
||||
0.d0, R(1,1), size(R, 1))
|
||||
|
||||
do k = 1, n_states_diag
|
||||
do i = 1, N
|
||||
U0(i,k) = R(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call ortho_qr(U0(1,1), size(U0, 1), N, n_states_diag)
|
||||
|
||||
do j = 1, n_states_diag
|
||||
k = 1
|
||||
do while((k < N) .and. (U0(k,j) == 0.d0))
|
||||
k = k+1
|
||||
enddo
|
||||
if(U0(k,j) * R(k,j) < 0.d0) then
|
||||
do i = 1, N
|
||||
W0(i,j) = -W0(i,j)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
enddo ! loop over while
|
||||
|
||||
! ---
|
||||
|
||||
write_buffer = '====='
|
||||
do i = 1, n_states
|
||||
write_buffer = trim(write_buffer)//' ================ ==========='
|
||||
enddo
|
||||
write(6,'(A)') trim(write_buffer)
|
||||
write(6,'(A)') ''
|
||||
|
||||
|
||||
print*, " Davidson eigenvalues"
|
||||
do k = 1, n_states
|
||||
Om(k) = h_val(k)
|
||||
print*, k, Om(k)
|
||||
enddo
|
||||
|
||||
deallocate(write_buffer)
|
||||
deallocate(i_omax)
|
||||
deallocate(to_print)
|
||||
deallocate(H_diag)
|
||||
deallocate(U0, U1)
|
||||
deallocate(W0, W1)
|
||||
deallocate(h)
|
||||
deallocate(h_vec)
|
||||
deallocate(h_val)
|
||||
deallocate(residual_norm)
|
||||
|
||||
if(kernel_name .eq. "gw") then
|
||||
deallocate(rho_tmp)
|
||||
deallocate(Om_tmp)
|
||||
endif
|
||||
|
||||
call wall_time(t2)
|
||||
write(*,'(A50, F12.4)') 'total wall time for Davidson (sec): ', t2-t1
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
@ -1,11 +1,12 @@
|
||||
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
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, &
|
||||
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_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)
|
||||
|
||||
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) :: doG0F2,doevGF2,doqsGF2
|
||||
logical,intent(in) :: doG0W0,doevGW,doqsGW
|
||||
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp
|
||||
|
||||
integer,intent(in) :: nNuc,nBas
|
||||
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
|
||||
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) :: doACFDT,exchange_kernel,doXBS
|
||||
|
||||
! 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_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_GF ,end_GF ,t_GF
|
||||
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 :: EGHF
|
||||
@ -290,4 +298,25 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
|
||||
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
|
||||
|
||||
@ -249,15 +249,15 @@ program QuAcK
|
||||
!--------------------------!
|
||||
! Generalized QuAcK branch !
|
||||
!--------------------------!
|
||||
|
||||
if(doGQuAcK) &
|
||||
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
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, &
|
||||
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_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)
|
||||
|
||||
!-----------!
|
||||
|
||||
Loading…
Reference in New Issue
Block a user