4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
This commit is contained in:
Pierre-Francois Loos 2021-10-19 18:07:44 +02:00
parent 87d4c8ac9a
commit 7fd40afee5
4 changed files with 220 additions and 3 deletions

View File

@ -13,9 +13,9 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
F F F F F
T F F F F
# G0T0 evGT qsGT
T F F
F F F
# MCMP2
F
# * unrestricted version available

View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000000001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
T T F T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0

215
src/MBPT/ufBSE.f90 Normal file
View File

@ -0,0 +1,215 @@
subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
! Unfold BSE@GW equations
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
! Local variables
integer :: s
integer :: i,j,k,l
integer :: a,b,c,d
integer :: ia,jb,kc,iajb,kcld
double precision :: tmp
integer :: n1h1p,n2h2p,nH
double precision,external :: Kronecker_delta
integer,allocatable :: order(:)
double precision,allocatable :: H(:,:)
double precision,allocatable :: X(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: Z(:)
! Output variables
! Hello world
write(*,*)
write(*,*)'**********************************************'
write(*,*)'| Unfolded BSE@GW calculation |'
write(*,*)'**********************************************'
write(*,*)
! TDA for W
write(*,*) 'Tamm-Dancoff approximation by default!'
write(*,*)
! Dimension of the supermatrix
n1h1p = nO*nV
n2h2p = nO*nO*nV*nV
nH = n1h1p + n2h2p + n2h2p
! Memory allocation
allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH))
! Initialization
H(:,:) = 0d0
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | A -Ve -Vh | !
! | | !
! H = | +Vh C2h2p 0 | !
! | | !
! | +Ve 0 C2p2h | !
! !
!---------------------------!
!---------!
! Block A !
!---------!
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ 2d0*ERI(i,b,a,j) - ERI(i,b,j,a)
end do
end do
end do
end do
!----------------!
! Blocks Vp & Ve !
!----------------!
iajb=0
do i=nC+1,nO
do a=nO+1,nBas-nR
do j=nC+1,nO
do b=nO+1,nBas-nR
iajb = iajb + 1
kc = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
kc = kc + 1
tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i)
H(n1h1p+iajb,kc ) = +tmp
H(kc ,n1h1p+n2h2p+iajb) = -tmp
tmp = sqrt(2d0)*Kronecker_delta(b,c)*ERI(a,k,i,j)
H(n1h1p+n2h2p+iajb,kc ) = +tmp
H(kc ,n1h1p+iajb) = -tmp
end do
end do
end do
end do
end do
end do
!-------------!
! Block 2h2p !
!-------------!
iajb = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do j=nC+1,nO
do b=nO+1,nBas-nR
iajb = iajb + 1
kcld = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do l=nC+1,nO
do d=nO+1,nBas-nR
kcld = kcld + 1
tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
+ 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d)
H(n1h1p +iajb,n1h1p +kcld) = tmp
H(n1h1p+n2h2p+iajb,n1h1p+n2h2p+kcld) = tmp
end do
end do
end do
end do
end do
end do
end do
end do
!-------------------------!
! Diagonalize supermatrix !
!-------------------------!
! call matout(nH,nH,H)
call diagonalize_general_matrix(nH,H,Om,X)
do s=1,nH
order(s) = s
end do
call quick_sort(Om,order,nH)
call set_order(X,order,nH,nH)
!-----------------!
! Compute weights !
!-----------------!
Z(:) = 0d0
do s=1,nH
do ia=1,n1h1p
Z(s) = Z(s) + X(ia,s)**2
end do
end do
!--------------!
! Dump results !
!--------------!
write(*,*)'-------------------------------------------'
write(*,*)' unfolded BSE excitation energies (eV) '
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','Omega (eV)','|','Z','|'
write(*,*)'-------------------------------------------'
do s=1,nH
if(Z(s) > 1d-7) &
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',Om(s)*HaToeV,'|',Z(s),'|'
enddo
write(*,*)'-------------------------------------------'
write(*,*)
end subroutine ufBSE

View File

@ -956,6 +956,8 @@ program QuAcK
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
! call ufBSE(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
end if
call cpu_time(end_G0W0)