PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 15:14:55 +02:00
Code Releases Activity

Added SRG for complex GW methods

Browse Source
This commit is contained in:
Loris Burth 2025-04-22 19:58:53 +02:00
parent 62f331af9e
commit 179a0835b3
10 changed files with 310 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/GW/RGW_self_energy.f90
View File

@ -149,5 +149,4 @@ subroutine RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
end do
end do
end do
end subroutine

9
src/GW/complex_RGW_QP_graph.f90
View File

@ -64,11 +64,16 @@ subroutine complex_RGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,Re_eHF,Im_eHF
nIt = nIt + 1
if(doSRG) then
call complex_RGW_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,nS,&
call complex_RGW_SRG_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,nS,&
Re_w,Im_w,Re_eOld,Im_eOld,Om,rho,&
Re_SigC,Im_SigC,Re_dSigC,Im_dSigC)
else
call complex_RGW_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,nS,&
Re_w,Im_w,Re_eOld,Im_eOld,Om,rho,&
Re_SigC,Im_SigC,Re_dSigC,Im_dSigC)
end if
Re_f = Re_w - Re_eHF(p) - Re_SigC
Im_f = Im_w - Im_eHF(p) - Im_SigC
Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)

88
src/GW/complex_RGW_SRG_SigC_dSigC.f90 Normal file
View File

@ -0,0 +1,88 @@
subroutine complex_RGW_SRG_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,nS,Re_w,Im_w,Re_e,Im_e,Om,rho,Re_SigC,Im_SigC,Re_DS,Im_DS)
! Complute diagonal of the correlation part of the self-energy and its derivative fully complex
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: eta
double precision,intent(in) :: flow
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: Re_e(nBas)
double precision,intent(in) :: Im_e(nBas)
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
complex*16,intent(in) :: Om(nS)
complex*16,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,m
double precision :: eps,s
double precision :: eta_tilde
complex*16 :: num
complex*16 :: tmp
! Output variables
double precision,intent(out) :: Re_SigC
double precision,intent(out) :: Im_SigC
double precision,intent(out) :: Re_DS
double precision,intent(out) :: Im_DS
! Initialize
Re_SigC = 0d0
Im_SigC = 0d0
Re_DS = 0d0
Im_DS = 0d0
s = flow
! Compute self energy and its derivative
! Occupied part
do i=nC+1,nO
do m=1,nS
eps = Re_w - Re_e(i) + real(Om(m))
eta_tilde = eta - Im_w + Im_e(i) - aimag(Om(m))
num = 2d0*rho(p,i,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC = Re_SigC + real(tmp)
Im_SigC = Im_SigC + aimag(tmp)
tmp = 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 = Re_DS + real(tmp)
Im_DS = Im_DS + aimag(tmp)
end do
end do
! Virtual part
do a=nO+1,nBas-nR
do m=1,nS
eps = Re_w - Re_e(a) - real(Om(m))
eta_tilde = eta + Im_w - Im_e(a) - aimag(Om(m))
num = 2d0*rho(p,a,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC = Re_SigC + real(tmp)
Im_SigC = Im_SigC + aimag(tmp)
tmp = 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 = Re_DS + real(tmp)
Im_DS = Im_DS + aimag(tmp)
end do
end do
end subroutine

193
src/GW/complex_RGW_SRG_self_energy.f90 Normal file
View File

@ -0,0 +1,193 @@
subroutine complex_RGW_SRG_self_energy(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
! Compute correlation part of the self-energy and the renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
double precision,intent(in) :: flow
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) :: e(nOrb)
complex*16,intent(in) :: Om(nS)
complex*16,intent(in) :: rho(nOrb,nOrb,nS)
! Local variables
integer :: i,j,a,b
integer :: p,q,m
double precision :: eps_p,eps_q,eta_tilde_p,eta_tilde_q
double precision :: eps,eta_tilde,s
complex*16 :: num,tmp
double precision,allocatable :: Re_DS(:)
double precision,allocatable :: Im_DS(:)
double precision,allocatable :: Re_Sig(:,:)
double precision,allocatable :: Im_Sig(:,:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
! Output variables
complex*16,intent(out) :: EcGM
complex*16,intent(out) :: Sig(nOrb,nOrb)
complex*16,intent(out) :: Z(nOrb)
!----------------!
! GW self-energy !
!----------------!
allocate(Re_DS(nOrb),Im_DS(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb))
Re_Sig(:,:) = 0d0
Im_Sig(:,:) = 0d0
Re_DS(:) = 0d0
Im_DS(:) = 0d0
s = flow
! Occupied part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) &
!$OMP PRIVATE(m,i,q,p,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,tmp) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
do m=1,nS
do i=nC+1,nO
eps_p = real(e(p)) - real(e(i)) + real(Om(m))
eps_q = real(e(q)) - real(e(i)) + real(Om(m))
eta_tilde_p = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m))
eta_tilde_q = eta - aimag(e(q)) + aimag(e(i)) - aimag(Om(m))
num = 2d0*rho(p,i,m)*rho(q,i,m)*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2)))
tmp = num*cmplx((eps_p + eps_q)/(eps_p**2 + eps_q**2 + eta_tilde_p**2 + eta_tilde_q**2),&
(eta_tilde_p + eta_tilde_q)/(eps_p**2 + eps_q**2 + eta_tilde_p**2 + eta_tilde_q**2),kind=8)
Re_Sig(p,q) = Re_Sig(p,q) + real(tmp)
Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp)
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(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) &
!$OMP PRIVATE(m,a,q,p,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,tmp) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nOrb-nR
do p=nC+1,nOrb-nR
do m=1,nS
do a=nO+1,nOrb-nR
eps_p = real(e(p)) - real(e(a)) - real(Om(m))
eps_q = real(e(q)) - real(e(a)) - real(Om(m))
eta_tilde_p = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m))
eta_tilde_q = eta + aimag(e(q)) - aimag(e(a)) - aimag(Om(m))
num = 2d0*rho(p,a,m)*rho(q,a,m)*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2)))
tmp = num*cmplx((eps_p +eps_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),&
-(eta_tilde_p + eta_tilde_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),kind=8)
Re_Sig(p,q) = Re_Sig(p,q) + real(tmp)
Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp)
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
!------------------------!
! Renormalization factor !
!------------------------!
!Occupied part of the renormalization factor
!$OMP PARALLEL &
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) &
!$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE)
!$OMP DO
do p=nC+1,nOrb-nR
do m=1,nS
do i=nC+1,nO
eps = real(e(p)) - real(e(i)) + real(Om(m))
eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m))
num = 2d0*rho(p,i,m)*rho(p,i,m)*(1d0-exp(-2d0*s*(eps**2+eta_tilde**2)))
tmp = 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(tmp)
Im_DS(p) = Im_DS(p) + aimag(tmp)
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the renormalization factor
!$OMP PARALLEL &
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) &
!$OMP PRIVATE(m,a,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE)
!$OMP DO
do p=nC+1,nOrb-nR
do m=1,nS
do a=nO+1,nOrb-nR
eps = real(e(p)) - real(e(a)) - real(Om(m))
eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m))
num = 2d0*rho(p,a,m)*rho(p,a,m)*(1d0-exp(-2d0*s*(eps**2+eta_tilde**2)))
tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
2*eta_tilde*eps/eps/(eps**2 + eta_tilde**2)**2,kind=8)
Re_DS(p) = Re_DS(p) + real(tmp)
Im_DS(p) = Im_DS(p) + aimag(tmp)
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Compute renormalization factor from derivative
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Z = cmplx(Re_Z,Im_Z,kind=8)
Sig = cmplx(Re_Sig,Im_Sig,kind=8)
deallocate(Re_DS)
deallocate(Im_DS)
deallocate(Re_Z)
deallocate(Im_Z)
deallocate(Re_Sig)
deallocate(Im_Sig)
!!-------------------------------------!
!! Galitskii-Migdal correlation energy !
!!-------------------------------------!
!
! EcGM = 0d0
! do m=1,nS
! do a=nO+1,nOrb-nR
! do i=nC+1,nO
!
! eps = e(a) - e(i) + Om(m)
! num = 4d0*rho(a,i,m)*rho(a,i,m)
! EcGM = EcGM - num*eps/(eps**2 + eta**2)
!
! end do
! end do
! end do
!
end subroutine

8
src/GW/complex_RGW_SRG_self_energy_diag.f90
View File

@ -24,7 +24,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re
! Local variables
integer :: i,a,p,m
double precision :: eps
double precision :: eps,s
complex*16 :: num
double precision :: eta_tilde
double precision,allocatable :: Re_DS(:)
@ -46,6 +46,8 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re
Re_DS(:) = 0d0
Im_DS(:) = 0d0
s = flow
!----------------!
! GW self-energy !
!----------------!
@ -56,7 +58,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re
do m=1,nS
eps = Re_e(p) - Re_e(i) + real(Om(m))
eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m))
num = 2d0*rho(p,i,m)**2
num = 2d0*rho(p,i,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
eta_tilde/(eps**2+eta_tilde**2),kind=8)
Re_Sig(p) = Re_Sig(p) + real(tmp)
@ -78,7 +80,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re
eps = Re_e(p) - Re_e(a) - real(Om(m))
eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m))
num = 2d0*rho(p,a,m)**2
num = 2d0*rho(p,a,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_Sig(p) = Re_Sig(p) + real(tmp)

2
src/GW/complex_RGW_self_energy.f90
View File

@ -41,7 +41,7 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Si
!----------------!
! GW self-energy !
!----------------!
allocate(Re_DS(nBas),Im_DS(nBas),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb))
allocate(Re_DS(nOrb),Im_DS(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb))
Re_Sig(:,:) = 0d0
Im_Sig(:,:) = 0d0

7
src/GW/complex_cRG0W0.f90
View File

@ -94,9 +94,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
flow = 500d0
if(doSRG) then
! Not implemented
write(*,*) '*** SRG regularized G0W0 scheme ***'
write(*,*) '!!! No SRG with cRG0W0 !!!'
write(*,*)
end if
@ -126,8 +124,11 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
!------------------------!
! Compute GW self-energy !
!------------------------!
if(doSRG) then
call complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
else
call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
end if
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!

9
src/GW/complex_evRGW.f90
View File

@ -146,10 +146,13 @@ subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
! Compute correlation part of the self-energy
! Implement here the srg version if(doSRG) .. complex_RGW_SRG_self_energy_diag
call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,&
if(doSRG) then
call complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,&
EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
else
call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,&
EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
end if
! Solve the quasi-particle equation

4
src/GW/complex_qsRGW.f90
View File

@ -215,8 +215,8 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
if(doSRG) then
write(*,*) "SRG not implemented"
!call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call complex_RGW_SRG_self_energy(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,&
rho,EcGM,SigC,Z)
else
call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,&
EcGM,SigC,Z)