10
1
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:24:58 +02:00

added complex G0F2

This commit is contained in:
Loris Burth 2025-03-25 15:36:04 +01:00
parent 738cdde624
commit 0ed01dc905
3 changed files with 299 additions and 0 deletions

153
src/GF/cRG0F2.f90 Normal file
View File

@ -0,0 +1,153 @@
subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: dophBSE
logical,intent(in) :: doppBSE
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nC
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) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: CAP(nOrb,nOrb)
! Local variables
integer :: p
double precision :: Ec
double precision :: EcBSE(nspin)
double precision,allocatable :: Re_SigC(:)
double precision,allocatable :: Im_SigC(:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
double precision,allocatable :: Re_eGFlin(:)
double precision, allocatable :: Im_eGFlin(:)
double precision,allocatable :: Re_eGF(:)
double precision,allocatable :: Im_eGF(:)
double precision, allocatable :: e_CAP(:)
! Hello world
write(*,*)
write(*,*)'*******************************'
write(*,*)'* Restricted G0F2 Calculation *'
write(*,*)'*******************************'
write(*,*)
! Memory allocation
allocate(Re_SigC(nOrb),Im_SigC(nOrb), Re_Z(nOrb),Im_Z(nOrb),&
Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),e_CAP(nOrb))
do p = 1, nOrb
e_CAP(p) = CAP(p,p)
end do
! Frequency-dependent second-order contribution
if(regularize) then
write(*,*) "Regularisation not implemented (yet)"
! call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
else
call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_CAP)
end if
Re_eGFlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGFlin(:) = e_CAP(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
if(linearize) then
write(*,*) '*** Quasiparticle energies obtained by linearization ***'
Re_eGF(:) = Re_eGFlin(:)
Im_eGF(:) = Im_eGFlin(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** NOT IMPLEMEMTED YET'
write(*,*)
!call RGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
end if
! Print results
! call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_cRG0F2(nOrb,nO,eHF,e_CAP,Re_SigC,Im_SigC,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec)
! Perform BSE@GF2 calculation
!
! if(dophBSE) then
!
! call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1)
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE)
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Perform ppBSE@GF2 calculation
!
! if(doppBSE) then
!
! call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
!
! EcBSE(2) = 3d0*EcBSE(2)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Testing zone
!
! if(dotest) then
!
! call dump_test_value('R','G0F2 correlation energy',Ec)
! call dump_test_value('R','G0F2 HOMO energy',eGF(nO))
! call dump_test_value('R','G0F2 LUMO energy',eGF(nO+1))
!
! end if
!
! deallocate(SigC, Z, eGFlin, eGF)
!
end subroutine

View File

@ -0,0 +1,88 @@
subroutine cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap)
! Compute diagonal part of the GF2 self-energy and its renormalization factor
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
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
integer :: p
double precision :: eps
double precision :: eta_tilde
double precision :: num
double precision,allocatable :: Re_DS(:)
double precision,allocatable :: Im_DS(:)
! Output variables
double precision,intent(out) :: Re_SigC(nBas)
double precision,intent(out) :: Im_SigC(nBas)
double precision,intent(out) :: Re_Z(nBas)
double precision,intent(out) :: Im_Z(nBas)
! Initialize
allocate(Re_DS(nBas),Im_DS(nBas))
Re_SigC(:) = 0d0
Im_SigC(:) = 0d0
Re_DS(:) = 0d0
Im_DS(:) = 0d0
! Compute GF2 self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = e(p) + e(a) - e(i) - e(j)
eta_tilde = eta - e_cap(p) + e_cap(i) + e_cap(a) - e_cap(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2)
Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(p) + e(i) - e(a) - e(b)
eta_tilde = eta + e_cap(p) - e_cap(a) - e_cap(b) + e_cap(i)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2)
Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2)
Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do
end do
end do
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
deallocate(Re_DS,Im_DS)
end subroutine

58
src/GF/print_cRG0F2.f90 Normal file
View File

@ -0,0 +1,58 @@
subroutine print_cRG0F2(nBas,nO,eHF,e_cap,Re_Sig,Im_Sig,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec)
! Print one-electron energies and other stuff for G0F2
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: e_cap(nBas)
double precision,intent(in) :: Re_Sig(nBas)
double precision,intent(in) :: Im_Sig(nBas)
double precision,intent(in) :: Re_eGF(nBas)
double precision,intent(in) :: Im_eGF(nBas)
double precision,intent(in) :: Re_Z(nBas)
double precision,intent(in) :: Im_Z(nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: Ec
integer :: p
integer :: HOMO
integer :: LUMO
double precision :: Gap
! Dump results
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)' One-shot G0F2 calculation'
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Re(Sig_GF2) (eV)','|','Re(Z)','|','Re(e_GF2) (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',Re_Sig(p)*HaToeV,'|',Re_Z(p),'|',Re_eGF(p)*HaToeV,'|'
end do
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','Im(e_HF) (eV)','|','Im(Sig_GF2) (eV)','|','Im(Z)','|','Im(e_GF2) (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',e_cap(p)*HaToeV,'|',Im_Sig(p)*HaToeV,'|',Im_Z(p),'|',Im_eGF(p)*HaToeV,'|'
end do
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A60,F15.6,A3)') 'G0F2 total energy =',ENuc + ERHF + Ec,' au'
write(*,'(2X,A60,F15.6,A3)') 'G0F2 correlation energy =',Ec,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine