4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 04:14:26 +01:00
This commit is contained in:
Pierre-Francois Loos 2024-09-11 17:21:43 +02:00
parent da7bdcd809
commit c2d92e2299
9 changed files with 233 additions and 109 deletions

View File

@ -24,7 +24,6 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
integer :: p,q
integer :: m
double precision :: Dpim,Dqim,Dpam,Dqam,Diam
double precision :: renorm
double precision :: s
! Output variables
@ -49,7 +48,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,i,q,p,Dpim,Dqim,renorm) &
!$OMP PRIVATE(m,i,q,p,Dpim,Dqim) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nOrb-nR
@ -59,8 +58,9 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
Dpim = e(p) - e(i) + Om(m)
Dqim = e(q) - e(i) + Om(m)
renorm = (1d0-exp(-s*Dpim*Dpim)*exp(-s*Dqim*Dqim))*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*renorm
SigC(p,q) = SigC(p,q) &
+ 2d0*rho(p,i,m)*rho(q,i,m)* &
(1d0-exp(-s*Dpim*Dpim)*exp(-s*Dqim*Dqim))*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
end do
end do
@ -73,7 +73,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) &
!$OMP PRIVATE(m,a,q,p,Dpam,Dqam,renorm) &
!$OMP PRIVATE(m,a,q,p,Dpam,Dqam) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nOrb-nR
@ -83,8 +83,9 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
Dpam = e(p) - e(a) - Om(m)
Dqam = e(q) - e(a) - Om(m)
renorm = (1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam))*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*renorm
SigC(p,q) = SigC(p,q) &
+ 2d0*rho(p,a,m)*rho(q,a,m)* &
(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam))*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
end do
end do

View File

@ -23,8 +23,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
integer :: i,j,a,b
integer :: p
integer :: m
double precision :: Dpim,Dqim,Dpam,Dqam,Diam
double precision :: renorm
double precision :: Dpim,Dpam,Diam
double precision :: s
! Output variables
@ -49,7 +48,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,i,p,Dpim,Dqim,renorm) &
!$OMP PRIVATE(m,i,p,Dpim) &
!$OMP DEFAULT(NONE)
!$OMP DO
do p=nC+1,nOrb-nR
@ -57,8 +56,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
do i=nC+1,nO
Dpim = e(p) - e(i) + Om(m)
renorm = (1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
SigC(p) = SigC(p) + 2d0*rho(p,i,m)**2*renorm
SigC(p) = SigC(p) + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
end do
end do
@ -70,7 +68,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) &
!$OMP PRIVATE(m,a,p,Dpam,Dqam,renorm) &
!$OMP PRIVATE(m,a,p,Dpam) &
!$OMP DEFAULT(NONE)
!$OMP DO
do p=nC+1,nOrb-nR
@ -78,8 +76,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
do a=nO+1,nOrb-nR
Dpam = e(p) - e(a) - Om(m)
renorm = (1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
SigC(p) = SigC(p) + 2d0*rho(p,a,m)**2*renorm
SigC(p) = SigC(p) + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
end do
end do

View File

@ -217,13 +217,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is))
end do
end if
call USRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z)
call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis

View File

@ -1,5 +1,5 @@
subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform unrestricted G0W0 calculation
@ -24,7 +24,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
@ -127,14 +127,12 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is),Om,rho(:,:,:,is))
end do
if(doSRG) then
call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
else
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
end if
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!

View File

@ -1,4 +1,4 @@
subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
! Compute correlation part of the self-energy
@ -21,10 +21,9 @@ subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
integer :: ispin
integer :: i,j,a,b
integer :: p,q,r
integer :: p,q
integer :: m
double precision :: Dpim,Dqim,Dpam,Dqam,Diam
double precision :: t1,t2
double precision :: s
! Output variables
@ -47,97 +46,95 @@ subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
! Occupied part of the correlation self-energy
! call wall_time(t1)
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(ispin,m,i,q,p,Dpim,Dqim) &
!$OMP DEFAULT(NONE)
!$OMP DO
do ispin=1,nspin
do q=nC(ispin)+1,nBas-nR(ispin)
do p=nC(ispin)+1,nBas-nR(ispin)
do q=nC(ispin)+1,nBas-nR(ispin)
do p=nC(ispin)+1,nBas-nR(ispin)
do m=1,nS
do i=nC(ispin)+1,nO(ispin)
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
Dqim = e(q,ispin) - e(i,ispin) + Om(m)
SigC(p,q,ispin) = SigC(p,q,ispin) &
+ rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) &
*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
end do
do i=nC(ispin)+1,nO(ispin)
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
Dqim = e(q,ispin) - e(i,ispin) + Om(m)
SigC(p,q,ispin) = SigC(p,q,ispin) &
+ rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) &
*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
end do
end do
end do
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! call wall_time(t2)
! print *, "first loop", (t2-t1)
! Virtual part of the correlation self-energy
! Virtual part of the correlation self-energy
call wall_time(t1)
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) &
!$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) &
!$OMP DEFAULT(NONE)
!$OMP DO
do ispin=1,nspin
do q=nC(ispin)+1,nBas-nR(ispin)
do p=nC(ispin)+1,nBas-nR(ispin)
do m=1,nS
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) &
!$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) &
!$OMP DEFAULT(NONE)
!$OMP DO
do ispin=1,nspin
do q=nC(ispin)+1,nBas-nR(ispin)
do p=nC(ispin)+1,nBas-nR(ispin)
do m=1,nS
do a=nO(ispin)+1,nBas-nR(ispin)
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
Dqam = e(q,ispin) - e(a,ispin) - Om(m)
SigC(p,q,ispin) = SigC(p,q,ispin) &
+ rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) &
*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
Dqam = e(q,ispin) - e(a,ispin) - Om(m)
SigC(p,q,ispin) = SigC(p,q,ispin) &
+ rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) &
*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
end do
end do
end do
end do
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! call wall_time(t2)
! print *, "second loop", (t2-t1)
! Initialize
!------------------------!
! Renormalization factor !
!------------------------!
Z(:,:) = 0d0
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do i=nC(ispin)+1,nO(ispin)
do m=1,nS
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2
end do
end do
end do
end do
! Virtual part of the correlation self-energy
! Occupied part of the renormalization factor
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do a=nO(ispin)+1,nBas-nR(ispin)
do p=nC(ispin)+1,nBas-nR(ispin)
do i=nC(ispin)+1,nO(ispin)
do m=1,nS
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2
end do
end do
end do
end do
end do
end do
! Compute renormalization factor from derivative of SigC
! Virtual part of the renormalization factor
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do a=nO(ispin)+1,nBas-nR(ispin)
do m=1,nS
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2
end do
end do
end do
end do
Z(:,:) = 1d0/(1d0 - Z(:,:))
! Galitskii-Migdal correlation energy
!-------------------------------------!
! Galitskii-Migdal correlation energy !
!-------------------------------------!
EcGM = 0d0
do ispin=1,nspin

View File

@ -0,0 +1,141 @@
subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
! Compute correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin)
! Local variables
integer :: ispin
integer :: i,j,a,b
integer :: p
integer :: m
double precision :: Dpim,Dpam,Diam
double precision :: s
! Output variables
double precision,intent(out) :: EcGM
double precision,intent(out) :: SigC(nBas,nspin)
double precision,intent(out) :: Z(nBas,nspin)
! SRG flow parameter
s = 500d0
! Initialize
SigC(:,:) = 0d0
!--------------------!
! SRG-GW self-energy !
!--------------------!
! Occupied part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(ispin,m,i,p,Dpim) &
!$OMP DEFAULT(NONE)
!$OMP DO
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do m=1,nS
do i=nC(ispin)+1,nO(ispin)
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
SigC(p,ispin) = SigC(p,ispin) + rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) &
!$OMP PRIVATE(ispin,m,a,p,Dpam) &
!$OMP DEFAULT(NONE)
!$OMP DO
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do m=1,nS
do a=nO(ispin)+1,nBas-nR(ispin)
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
SigC(p,ispin) = SigC(p,ispin) + rho(p,a,m,ispin)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
!------------------------!
! Renormalization factor !
!------------------------!
Z(:,:) = 0d0
! Occupied part of the renormalization factor
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do i=nC(ispin)+1,nO(ispin)
do m=1,nS
Dpim = e(p,ispin) - e(i,ispin) + Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2
end do
end do
end do
end do
! Virtual part of the renormalization factor
do ispin=1,nspin
do p=nC(ispin)+1,nBas-nR(ispin)
do a=nO(ispin)+1,nBas-nR(ispin)
do m=1,nS
Dpam = e(p,ispin) - e(a,ispin) - Om(m)
Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2
end do
end do
end do
end do
Z(:,:) = 1d0/(1d0 - Z(:,:))
!-------------------------------------!
! Galitskii-Migdal correlation energy !
!-------------------------------------!
EcGM = 0d0
do ispin=1,nspin
do i=nC(ispin)+1,nO(ispin)
do a=nO(ispin)+1,nBas-nR(ispin)
do m=1,nS
Diam = e(a,ispin) - e(i,ispin) + Om(m)
EcGM = EcGM - rho(a,i,m,ispin)*rho(a,i,m,ispin)*(1d0-exp(-2d0*s*Diam*Diam))/Diam
end do
end do
end do
end do
end subroutine

View File

@ -125,8 +125,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Compute correlation part of the self-energy
call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho)
if(doSRG) then
call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
else

View File

@ -1,5 +1,5 @@
subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA, &
spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, &
spin_conserved,spin_flip,linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform self-consistent eigenvalue-only GW calculation
@ -28,7 +28,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
@ -153,13 +153,12 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is))
end do
if(doSRG) then
call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
else
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
end if
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
!-----------------------------------!
! Solve the quasi-particle equation !

View File

@ -1,5 +1,5 @@
subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -25,7 +25,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -217,13 +217,12 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is))
end do
if(doSRG) then
call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
else
call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
end if
call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
! Make correlation self-energy Hermitian and transform it back to AO basis