4
1
mirror of https://github.com/pfloos/quack synced 2025-05-07 07:44:43 +02:00

start cRG0W0

This commit is contained in:
Loris Burth 2025-03-14 09:50:54 +01:00
parent 9498174b75
commit 0b8e159b28
6 changed files with 121 additions and 10 deletions

View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF HFB cRHF
F F F F F T
T F F F F F
# MP2 MP3
F F
# CCD pCCD DCD CCSD CCSD(T)
@ -12,11 +12,13 @@
F F F F
# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
F F F F F F
# G0W0 evGW qsGW ufG0W0 ufGW
F F F F F
# G0W0 evGW qsGW ufG0W0 ufGW
F F F F F
# G0T0pp evGTpp qsGTpp ufG0T0pp
F F F F
# G0T0eh evGTeh qsGTeh
F F F
# cG0W0
T
# Rtest Utest Gtest
F F F

View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS guess mix shift stab search
256 0.00001 5 1 0.0 0.0 F F
1000 0.00001 5 1 0.0 0.0 F F
# MP: reg
F
# CC: maxSCF thresh DIIS
@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg
256 0.00001 5 F 0.0 F F
256 0.00001 5 T 0.0 F F
# GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS

View File

@ -1,4 +1,4 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT, &
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
@ -17,6 +17,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
logical,intent(in) :: doqsGW
logical,intent(in) :: doufG0W0
logical,intent(in) :: doufGW
logical,intent(in) :: docG0W0
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -87,6 +88,22 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
end if
!------------------------------------------------------------------------
! Perform cG0W0 calculation
!------------------------------------------------------------------------
if(docG0W0) then
call wall_time(start_GW)
call cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGW calculation
!------------------------------------------------------------------------

View File

@ -1,4 +1,4 @@
subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform G0W0 calculation
@ -120,10 +120,12 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Compute GW self-energy !
!------------------------!
!!! HERE I HAVE TO WORK !!!
if(doSRG) then
call RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
write(*,*) "SRG for this method is not yet implemented !"
!call RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
else
call RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
call cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
end if
!-----------------------------------!

View File

@ -0,0 +1,91 @@
subroutine RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
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
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,p,m
double precision :: num,eps
! Output variables
double precision,intent(out) :: Sig(nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: EcGM
! Initialize
Sig(:) = 0d0
Z(:) = 0d0
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do m=1,nS
eps = e(p) - e(i) + Om(m)
num = 2d0*rho(p,i,m)**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
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do m=1,nS
eps = e(p) - e(a) - Om(m)
num = 2d0*rho(p,a,m)**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
end do
end do
! Galitskii-Migdal correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do m=1,nS
eps = e(a) - e(i) + Om(m)
num = 4d0*rho(a,i,m)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end do
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine

View File

@ -22,7 +22,6 @@ subroutine read_CAP_integrals(nBas,W)
! Open file with integrals
debug = .false.
open(unit=31,file='int/CAP.dat')
! Read CAP integrals