4
1
mirror of https://github.com/pfloos/quack synced 2024-12-31 08:36:05 +01:00
This commit is contained in:
Pierre-Francois Loos 2024-10-30 10:08:39 +01:00
parent df7f38b7b5
commit 1a4ad3e225
7 changed files with 55 additions and 53 deletions

View File

@ -176,13 +176,12 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
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)
call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipole_int,eHF,eGT,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy = ',EcBSE,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF total energy = ',ENuc + EGHF + EcBSE,' au'
write(*,'(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(*,*)

View File

@ -46,21 +46,26 @@ subroutine GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2
do r=nC+1,nOrb-nR
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r)
cd = 0
do c = nO+1, nOrb-nR
do d = c+1, nOrb-nR
do c=nO+1,nOrb-nR
do d=c+1,nOrb-nR
cd = cd + 1
T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd) * rho1(r,s,cd) / Om1(cd)
T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd)*rho1(r,s,cd)/Om1(cd)
end do ! d
end do ! c
kl = 0
do k = nC+1, nO
do l = k+1, nO
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl) * rho2(r,s,kl) / Om2(kl)
T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl)*rho2(r,s,kl)/Om2(kl)
enddo
enddo

View File

@ -1,5 +1,4 @@
subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
ERI,dipole_int,eT,eGT,EcBSE)
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
@ -30,7 +29,7 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
! Local variables
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
double precision,allocatable :: Om1(:), Om2(:)
double precision,allocatable :: X1(:,:), X2(:,:)
@ -46,14 +45,13 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
!----------------------------------------------
! Compute linear response
!----------------------------------------------
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
Bpp(:,:) = 0d0
Cpp(:,:) = 0d0
Dpp(:,:) = 0d0
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eT,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eT,ERI,Dpp)
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppGLR_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)
@ -62,7 +60,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
!----------------------------------------------
! 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)
@ -73,9 +73,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp)
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
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
@ -91,12 +91,10 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KC_sta)
call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KD_sta)
if(.not.TDA_T) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, &
Om1,rho1,Om2,rho2,T,KB_sta)
call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KC_sta)
call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KD_sta)
if(.not.TDA_T) &
call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta)
deallocate(Om1,Om2,rho1,rho2)
! Deallocate the 4-tensor T

View File

@ -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

View File

@ -48,14 +48,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1
do cd=1,nVVs
eps = e(p) + e(i) - Om1s(cd)
num = (1d0/2d0)*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 = (3d0/2d0)*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 = (1d0/2d0)*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 = (3d0/2d0)*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 = (1d0/2d0)*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 = (3d0/2d0)*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 = (1d0/2d0)*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 = (3d0/2d0)*rho2t(a,b,kl)**2
num = 1.5d0*rho2t(a,b,kl)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do

View File

@ -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(*,*)

View File

@ -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(*,*)