10
1
mirror of https://github.com/pfloos/quack synced 2025-05-07 07:35:02 +02:00

first buggy version of complex G0F2

This commit is contained in:
Loris Burth 2025-03-27 15:35:13 +01:00
parent 35c2864ea7
commit 8ceaecd4d7
7 changed files with 349 additions and 14 deletions

View File

@ -90,7 +90,7 @@ subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** NOT IMPLEMEMTED YET'
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,eHF,e_cap,Re_eGF,Im_eGF,Re_Z,Im_Z)

View File

@ -56,7 +56,7 @@ subroutine cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Re_SigC,Im_SigC,Re_
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
Im_DS(p) = Im_DS(p) - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
end do
end do

81
src/GF/complex_RGF.f90 Normal file
View File

@ -0,0 +1,81 @@
subroutine complex_RGF(dotest,docG0F2,maxSCF, &
thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a one-shot second-order Green function calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: docG0F2
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
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) :: 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) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
complex*16,intent(in) :: ERHF
complex*16,intent(in) :: eHF(nOrb)
complex*16,intent(in) :: cHF(nBas,nOrb)
complex*16,intent(in) :: PHF(nBas,nBas)
complex*16,intent(in) :: S(nBas,nBas)
complex*16,intent(in) :: CAP(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
complex*16,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
complex*16,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
double precision :: start_GF ,end_GF ,t_GF
!------------------------------------------------------------------------
! Compute complex G0F2 electronic binding energies
!------------------------------------------------------------------------
if(docG0F2) then
call wall_time(start_GF)
call complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
end subroutine

98
src/GF/complex_cRG0F2.f90 Normal file
View File

@ -0,0 +1,98 @@
subroutine complex_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
complex*16,intent(in) :: eHF(nOrb)
complex*16,intent(in) :: CAP(nOrb,nOrb)
complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
complex*16,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! 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 :: Re_eHF(:)
double precision, allocatable :: Im_eHF(:)
! 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),Re_eHF(nOrb),Im_eHF(nOrb))
Re_eHF(:) = real(eHF(:))
Im_eHF(:) = aimag(eHF(:))
! Frequency-dependent second-order contribution
call cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,real(eHF),aimag(eHF),ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
Re_eGFlin(:) = Re_eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGFlin(:) = Im_eHF(:) + 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 *** '
write(*,*)
write(*,*) "ONLY LINEARISATION IMPLEMENTED YET"
!call cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,eHF,e_cap,Re_eGF,Im_eGF,Re_Z,Im_Z)
end if
! Print results
call print_complex_cRG0F2(nOrb,nO,Re_eHF,Im_eHF,Re_SigC,Im_SigC,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec)
deallocate(Re_SigC,Im_SigC, Re_Z,Im_Z,&
Re_eGFlin,Im_eGFlin, Re_eGF,Im_eGF)
end subroutine

View File

@ -0,0 +1,97 @@
subroutine complex_cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,Re_e,Im_e,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
! 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) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
integer :: p
double precision :: eps
double precision :: eta_tilde
complex*16 :: num
double precision,allocatable :: Re_DS(:)
double precision,allocatable :: Im_DS(:)
complex*16 :: z_dummy
! 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 = Re_e(p) + Re_e(a) - Re_e(i) - Re_e(j)
eta_tilde = eta - Im_e(p) + Im_e(i) + Im_e(a) - Im_e(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC(p) = Re_SigC(p) + real(z_dummy)
Im_SigC(p) = Im_SigC(p) + aimag(z_dummy)
write(*,*) "Numerator"
write(*,*) num
z_dummy = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
-2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8)
Re_DS(p) = Re_DS(p) + real(z_dummy)
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
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 = Re_e(p) + Re_e(i) - Re_e(a) - Re_e(b)
eta_tilde = eta + Im_e(p) - Im_e(a) - Im_e(b) + Im_e(i)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC(p) = Re_SigC(p) + real(z_dummy)
Im_SigC(p) = Im_SigC(p) + aimag(z_dummy)
write(*,*) "Numerator"
write(*,*) num
z_dummy = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8)
Re_DS(p) = Re_DS(p) + real(z_dummy)
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
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

View File

@ -0,0 +1,47 @@
subroutine print_complex_cRG0F2(nBas,nO,Re_eHF,Im_eHF,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) :: Re_eHF(nBas)
double precision,intent(in) :: Im_eHF(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)') &
'|','#','|','Re(e_HF) (eV)','|','Re(Sig_GF2) (eV)','|','Re(Z)','|','Re(e_GF2) (eV)','|'
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,'|',Re_eHF(p)*HaToeV,'|',Re_Sig(p)*HaToeV,'|',Re_Z(p),'|',Re_eGF(p)*HaToeV,'|'
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,'|',Im_eHF(p)*HaToeV,'|',Im_Sig(p)*HaToeV,'|',Im_Z(p),'|',Im_eGF(p)*HaToeV,'|'
write(*,*)'-------------------------------------------------------------------------------'
end do
end subroutine

View File

@ -224,7 +224,6 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
if (doCAP) then
call complex_AOtoMO(nBas,nOrb,complex_cHF,CAP_AO,complex_CAP_MO)
complex_CAP_MO = (0d0,1d0)*complex_CAP_MO
call complex_matout(nBas,nOrb,complex_CAP_MO)
end if
else
@ -241,11 +240,6 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
end if
call wall_time(end_AOtoMO)
! Transform CAP integrals
if (doCAP .and. (.not. docRHF)) then
write(*,*) "blub 1"
end if
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
@ -365,11 +359,11 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 .or. docG0F2
if(doGF) then
if(doGF .and. .not. docRHF) then
call wall_time(start_GF)
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, &
thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GF)
@ -379,6 +373,24 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
end if
!---------------------------------!
! complex Green's function module !
!---------------------------------!
if(doGF .and. docRHF) then
call wall_time(start_GF)
call complex_RGF(dotest,docG0F2,maxSCF_GF, &
thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,complex_ERHF, &
S,X,T,V,Hc,ERI_AO,complex_ERI_MO,complex_CAP_MO,dipole_int_AO,complex_dipole_int_MO,&
complex_PHF,complex_cHF,complex_eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
!-----------!
! GW module !
!-----------!
@ -400,9 +412,9 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
end if
!------------!
! cGW module !
!------------!
!-------------------!
! complex GW module !
!-------------------!
! IMPLEMENT LATER TO TREAT EVERYTHING COMPLEX, i.e. start from complex HF orbitals