4
1
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:53 +01:00
quack/src/GT/ufRG0T0pp.f90

610 lines
16 KiB
Fortran
Raw Normal View History

2024-09-03 16:44:57 +02:00
subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
2023-11-28 19:15:39 +01:00
! Upfolded G0T0pp equations
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
2023-11-29 16:20:53 +01:00
logical,intent(in) :: TDA_T
2023-11-28 19:15:39 +01:00
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)
! Local variables
integer :: p
integer :: s
integer :: i,j,k,l
integer :: a,b,c,d
2023-11-29 16:20:53 +01:00
integer :: ij,ab
2023-11-28 19:15:39 +01:00
integer :: klc,kcd,ija,ijb,iab,jab
2023-11-29 17:16:33 +01:00
logical :: print_T = .false.
2023-11-28 19:15:39 +01:00
logical :: dRPA
integer :: ispin
2023-11-29 16:20:53 +01:00
integer :: iblock
2023-11-29 22:57:50 +01:00
integer :: nOOs,nOOt
integer :: nVVs,nVVt
2023-11-29 16:20:53 +01:00
double precision :: EcRPA(nspin)
2023-11-28 19:15:39 +01:00
integer :: n2h1p,n2p1h,nH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
2023-11-29 16:20:53 +01:00
double precision,allocatable :: eGT(:)
2023-11-28 19:15:39 +01:00
double precision,allocatable :: Z(:)
2023-11-29 16:20:53 +01:00
double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Cpp(:,:)
double precision,allocatable :: Dpp(:,:)
double precision,allocatable :: Om1s(:),Om1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Om2s(:),Om2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
2023-11-28 19:15:39 +01:00
2023-12-03 22:12:41 +01:00
logical :: verbose = .true.
2023-11-28 19:15:39 +01:00
double precision,parameter :: cutoff1 = 0.01d0
double precision,parameter :: cutoff2 = 0.01d0
double precision :: eF
2023-12-14 21:44:17 +01:00
double precision,parameter :: window = 1.5d0
2023-11-28 19:15:39 +01:00
double precision :: start_timing,end_timing,timing
! Output variables
! Hello world
write(*,*)
write(*,*)'******************************************'
write(*,*)'* Restricted Upfolded G0T0pp Calculation *'
write(*,*)'******************************************'
write(*,*)
2023-11-29 16:20:53 +01:00
! Dimensions of the ppRPA linear reponse matrices
2023-11-28 19:15:39 +01:00
2023-11-30 21:21:05 +01:00
! nOOs = nO*(nO + 1)/2
! nVVs = nV*(nV + 1)/2
2023-11-28 19:15:39 +01:00
2023-11-30 21:21:05 +01:00
nOOs = nO*nO
nVVs = nV*nV
2023-11-28 19:15:39 +01:00
2023-11-29 16:20:53 +01:00
nOOt = nO*(nO - 1)/2
nVVt = nV*(nV - 1)/2
2023-11-29 22:57:50 +01:00
! nOO = nO*nO
! nVV = nV*nV
2023-11-29 16:20:53 +01:00
! Dimension of the supermatrix
2023-11-29 17:21:36 +01:00
n2h1p = (nOOs+nOOt)*nV
n2p1h = (nVVs+nVVt)*nO
2023-11-29 16:20:53 +01:00
nH = 1 + n2h1p + n2p1h
2023-11-28 19:15:39 +01:00
! Initialization
dRPA = .true.
EcRPA = 0d0
eF = 0.5d0*(eHF(nO+1) + eHF(nO))
!------------------!
! Compute T-matrix !
!------------------!
if(.not. TDA_T) then
! Memory allocation
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
! alpha-beta block
ispin = 1
2023-11-30 21:21:05 +01:00
! iblock = 1
iblock = 3
! Compute linear response
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
2023-11-30 21:21:05 +01:00
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(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s(:))
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s(:))
! Compute excitation densities
2024-09-03 16:44:57 +02:00
call RGTpp_excitation_density(iblock,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
2023-11-30 21:21:05 +01:00
! iblock = 2
iblock = 4
! Compute linear response
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
2023-11-30 21:21:05 +01:00
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(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOt,Om2t)
! Compute excitation densities
2024-09-03 16:44:57 +02:00
call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
deallocate(Bpp,Cpp,Dpp,X1t,Y1t,X2t,Y2t)
end if
2023-11-29 16:20:53 +01:00
! Memory allocation
allocate(H(nH,nH),eGT(nH),Z(nH))
2023-11-28 19:15:39 +01:00
!-------------------------!
! Main loop over orbitals !
!-------------------------!
2023-11-29 22:57:50 +01:00
do p=nO,nO+1
2023-11-28 19:15:39 +01:00
H(:,:) = 0d0
2023-11-29 16:20:53 +01:00
if (TDA_T) then
2023-11-28 19:15:39 +01:00
2023-11-29 16:20:53 +01:00
! TDA for T
2023-11-28 19:15:39 +01:00
write(*,*) 'Tamm-Dancoff approximation actived!'
write(*,*)
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F V2h1p V2p1h | !
! | | !
! H = | V2h1p C2h1p 0 | !
! | | !
! | V2p1h 0 C2p1h | !
! !
!---------------------------!
call wall_time(start_timing)
!---------!
! Block F !
!---------!
! H(1,1) = eHF(p)
!-------------!
! Block V2h1p !
!-------------!
2023-11-28 19:15:39 +01:00
! ija = 0
! do i=nC+1,nO
! do j=nC+1,nO
! do a=nO+1,nBas-nR
! ija = ija + 1
!
! H(1 ,1+ija) = ERI(p,a,i,j)
! H(1+ija,1 ) = ERI(p,a,i,j)
2023-11-28 19:15:39 +01:00
!
! end do
! end do
! end do
! !-------------!
! ! Block V2p1h !
! !-------------!
!
! iab = 0
! do i=nC+1,nO
! do a=nO+1,nBas-nR
! do b=nO+1,nBas-nR
! iab = iab + 1
!
! H(1 ,1+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a)
! H(1+n2h1p+iab,1 ) = sqrt(2d0)*ERI(p,i,b,a)
!
! end do
! end do
! end do
!
! !-------------!
! ! Block C2h1p !
! !-------------!
!
! ija = 0
! do i=nC+1,nO
! do j=nC+1,nO
! do a=nO+1,nBas-nR
! ija = ija + 1
!
! klc = 0
! do k=nC+1,nO
! do l=nC+1,nO
! do c=nO+1,nBas-nR
! klc = klc + 1
!
! H(1+ija,1+klc) &
! = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
! - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k)
!
! end do
! end do
! end do
!
! end do
! end do
! end do
!
! !-------------!
! ! Block C2p1h !
! !-------------!
!
! iab = 0
! do i=nC+1,nO
! do a=nO+1,nBas-nR
! do b=nO+1,nBas-nR
! iab = iab + 1
!
! kcd = 0
! do k=nC+1,nO
! do c=nO+1,nBas-nR
! do d=nO+1,nBas-nR
! kcd = kcd + 1
!
! H(1+n2h1p+iab,1+n2h1p+kcd) &
! = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
! + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d)
!
! end do
! end do
! end do
!
! end do
! end do
! end do
!
call wall_time(end_timing)
timing = end_timing - start_timing
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds'
write(*,*)
else
2023-11-29 16:20:53 +01:00
! RPA for T
2023-11-28 19:15:39 +01:00
write(*,*) 'Tamm-Dancoff approximation deactivated!'
write(*,*)
!---------------------------!
! Compute GW supermatrix !
!---------------------------!
! !
! | F W2h1p W2p1h | !
! | | !
! H = | W2h1p D2h1p 0 | !
! | | !
! | W2p1h 0 D2p1h | !
! !
!---------------------------!
call wall_time(start_timing)
!---------!
! Block F !
!---------!
H(1,1) = eHF(p)
!-------------!
! Block D2h1p !
!-------------!
ija = 0
2023-11-29 16:20:53 +01:00
do ij=1,nOOs
do a=nO+1,nBas-nR
ija = ija + 1
2023-11-29 17:16:33 +01:00
H(1+ija,1+ija) = - eHF(a) + Om2s(ij)
2023-11-29 16:20:53 +01:00
end do
end do
do ij=1,nOOt
do a=nO+1,nBas-nR
2023-11-28 19:15:39 +01:00
ija = ija + 1
2023-11-29 17:16:33 +01:00
H(1+ija,1+ija) = - eHF(a) + Om2t(ij)
2023-11-28 19:15:39 +01:00
end do
end do
!-------------!
! Block W2h1p !
!-------------!
ija = 0
2023-11-29 16:20:53 +01:00
do ij=1,nOOs
do a=nO+1,nBas-nR
2023-11-28 19:15:39 +01:00
ija = ija + 1
2023-11-30 21:21:05 +01:00
H(1 ,1+ija) = rho2s(p,a,ij)
H(1+ija,1 ) = rho2s(p,a,ij)
2023-11-29 16:20:53 +01:00
end do
end do
do ij=1,nOOt
do a=nO+1,nBas-nR
ija = ija + 1
2023-11-30 21:21:05 +01:00
H(1 ,1+ija) = rho2t(p,a,ij)
H(1+ija,1 ) = rho2t(p,a,ij)
2023-11-28 19:15:39 +01:00
end do
end do
!-------------!
! Block D2p1h !
!-------------!
iab = 0
2023-11-29 16:20:53 +01:00
do ab=1,nVVs
do i=nC+1,nO
iab = iab + 1
H(1+n2h1p+iab,1+n2h1p+iab) = - eHF(i) + Om1s(ab)
end do
end do
do ab=1,nVVt
do i=nC+1,nO
2023-11-28 19:15:39 +01:00
iab = iab + 1
2023-11-29 17:16:33 +01:00
H(1+n2h1p+iab,1+n2h1p+iab) = - eHF(i) + Om1t(ab)
2023-11-28 19:15:39 +01:00
end do
end do
!-------------!
! Block W2p1h !
!-------------!
iab = 0
2023-11-29 16:20:53 +01:00
do ab=1,nVVs
do i=nC+1,nO
iab = iab + 1
2023-11-30 21:21:05 +01:00
H(1 ,1+n2h1p+iab) = rho1s(p,i,ab)
H(1+n2h1p+iab,1 ) = rho1s(p,i,ab)
2023-11-29 16:20:53 +01:00
end do
end do
do ab=1,nVVt
do i=nC+1,nO
2023-11-28 19:15:39 +01:00
iab = iab + 1
2023-11-30 21:21:05 +01:00
H(1 ,1+n2h1p+iab) = rho1t(p,i,ab)
H(1+n2h1p+iab,1 ) = rho1t(p,i,ab)
2023-11-28 19:15:39 +01:00
end do
end do
call wall_time(end_timing)
timing = end_timing - start_timing
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds'
write(*,*)
end if
!-------------------------!
! Diagonalize supermatrix !
!-------------------------!
call wall_time(start_timing)
2023-11-29 16:20:53 +01:00
call diagonalize_matrix(nH,H,eGT)
2023-11-28 19:15:39 +01:00
call wall_time(end_timing)
timing = end_timing - start_timing
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for diagonalization of supermatrix = ',timing,' seconds'
write(*,*)
!-----------------!
! Compute weights !
!-----------------!
do s=1,nH
Z(s) = H(1,s)**2
end do
!--------------!
! Dump results !
!--------------!
write(*,*)'-------------------------------------------'
2023-11-29 17:16:33 +01:00
write(*,'(1X,A34,I3,A6)')'| G0T0pp energies (eV) for orbital',p,' |'
2023-11-28 19:15:39 +01:00
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP','|','Z','|'
write(*,*)'-------------------------------------------'
do s=1,nH
2023-12-14 21:44:17 +01:00
if(eGT(s) < eF .and. eGT(s) > eF - window) then
!if(Z(s) > cutoff1) then
2023-11-28 19:15:39 +01:00
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
2023-11-29 16:20:53 +01:00
'|',s,'|',eGT(s)*HaToeV,'|',Z(s),'|'
2023-11-28 19:15:39 +01:00
end if
end do
write(*,*)'-------------------------------------------'
write(*,*)
if(verbose) then
2023-11-29 16:20:53 +01:00
if(TDA_T) then
2023-11-28 19:15:39 +01:00
! TDA printing format
2023-11-29 16:20:53 +01:00
! do s=1,nH
!
! if(eGT(s) < eF .and. eGT(s) > eF - window) then
!
! write(*,*)'-------------------------------------------------------------'
! write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') &
! 'Orbital',p,' and #',s,':','e_QP = ',eGT(s)*HaToeV,' eV and Z = ',Z(s)
! write(*,*)'-------------------------------------------------------------'
! write(*,'(1X,A20,1X,A20,1X,A15,1X)') &
! ' Configuration ',' Coefficient ',' Weight '
! write(*,*)'-------------------------------------------------------------'
!
! if(p <= nO) &
! write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
! ' (',p,') ',H(1,s),H(1,s)**2
! if(p > nO) &
! write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') &
! ' (',p,') ',H(1,s),H(1,s)**2
!
! ija = 0
! do i=nC+1,nO
! do j=nC+1,nO
! do a=nO+1,nBas-nR
! ija = ija + 1
!
! if(abs(H(1+ija,s)) > cutoff2) &
! write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') &
! ' (',i,',',j,') -> (',a,') ',H(1+ija,s),H(1+ija,s)**2
!
! end do
! end do
! end do
!
! iab = 0
! do i=nC+1,nO
! do a=nO+1,nBas-nR
! do b=nO+1,nBas-nR
! iab = iab + 1
!
! if(abs(H(1+n2h1p+iab,s)) > cutoff2) &
! write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') &
! ' (',i,') -> (',a,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2
!
! end do
! end do
! end do
2023-11-28 19:15:39 +01:00
2023-11-29 16:20:53 +01:00
! write(*,*)'-------------------------------------------------------------'
! write(*,*)
2023-11-28 19:15:39 +01:00
2023-11-29 16:20:53 +01:00
! end if
2023-11-28 19:15:39 +01:00
2023-11-29 16:20:53 +01:00
! end do
2023-11-28 19:15:39 +01:00
else
! non-TDA printing format
do s=1,nH
2023-12-14 11:16:44 +01:00
if(eGT(s) < eF .and. eGT(s) > eF - window) then
2023-12-14 21:44:17 +01:00
!if(Z(s) > cutoff2) then
2023-11-28 19:15:39 +01:00
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') &
2023-11-29 16:20:53 +01:00
'Orbital',p,' and #',s,':','e_QP = ',eGT(s)*HaToeV,' eV and Z = ',Z(s)
2023-11-28 19:15:39 +01:00
write(*,*)'-------------------------------------------------------------'
2023-11-29 17:16:33 +01:00
write(*,'(1X,A24,1X,A20,1X,A15,1X)') &
2023-11-29 16:20:53 +01:00
' Conf. (i,ab) or (a,ij) ',' Coefficient ',' Weight '
2023-11-28 19:15:39 +01:00
write(*,*)'-------------------------------------------------------------'
if(p <= nO) &
2023-12-15 13:50:05 +01:00
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',p,') ',H(1,s),H(1,s)**2,-eHF(p)*HaToeV
2023-11-28 19:15:39 +01:00
if(p > nO) &
2023-12-15 13:50:05 +01:00
write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',p,') ',H(1,s),H(1,s)**2,-eHF(p)*HaToeV
2023-11-28 19:15:39 +01:00
ija = 0
2023-11-29 22:57:50 +01:00
do ij=1,nOOs+nOOt
2023-11-29 16:20:53 +01:00
do a=nO+1,nBas-nR
2023-11-28 19:15:39 +01:00
ija = ija + 1
2023-12-15 13:50:05 +01:00
if(abs(H(1+ija,s)) > cutoff2 .and. ij<nOOs+1) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',a,',',ij,') ',H(1+ija,s),H(1+ija,s)**2,(- eHF(a) + Om2s(ij))*HaToeV
if(abs(H(1+ija,s)) > cutoff2 .and. ij>nOOs) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',a,',',ij,') ',H(1+ija,s),H(1+ija,s)**2,(- eHF(a) + Om2t(ij-nOOs))*HaToeV
2023-11-28 19:15:39 +01:00
end do
end do
iab = 0
2023-11-29 22:57:50 +01:00
do ab=1,nVVs+nVVt
2023-11-29 16:20:53 +01:00
do i=nC+1,nO
2023-11-28 19:15:39 +01:00
iab = iab + 1
2023-12-15 13:50:05 +01:00
if(abs(H(1+n2h1p+iab,s)) > cutoff2 .and. ab<nVVs+1) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',i,',',ab,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2,(-eHF(i) + Om1s(ab))*HaToeV
if(abs(H(1+n2h1p+iab,s)) > cutoff2 .and. ab>nVVs) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6,1X,F12.6)') &
' (',i,',',ab,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2,(-eHF(i) + Om1t(ab-nVVs))*HaToeV
2023-11-28 19:15:39 +01:00
end do
end do
write(*,*)'-------------------------------------------------------------'
write(*,*)
end if
end do
end if
end if
end do
end subroutine