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

regularization for complex GF2 self energy diag

Browse Source
This commit is contained in:
Loris Burth 2025-04-24 18:06:32 +02:00
parent 9a39e67a3b
commit cb647905cc
9 changed files with 213 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/GF/RGF2_reg_self_energy_diag.f90
View File

@ -81,5 +81,4 @@ subroutine RGF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
end do end do
Z(:) = 1d0/(1d0 - Z(:)) Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

1
src/GF/RGF2_self_energy_diag.f90
View File

@ -66,7 +66,6 @@ subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
end do end do
end do end do
end do end do
Z(:) = 1d0/(1d0 - Z(:)) Z(:) = 1d0/(1d0 - Z(:))
end subroutine end subroutine

3
src/GF/cRG0F2.f90
View File

@ -51,7 +51,6 @@ subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize
double precision, allocatable :: e_CAP(:) double precision, allocatable :: e_CAP(:)
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'*******************************' write(*,*)'*******************************'
write(*,*)'* Restricted G0F2 Calculation *' write(*,*)'* Restricted G0F2 Calculation *'
@ -70,7 +69,7 @@ subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize
if(regularize) then if(regularize) then
write(*,*) "Regularisation not implemented (yet)" write(*,*) "Regularisation not implemented (yet)"
! call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z) !call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
else else

12
src/GF/complex_cRG0F2.f90
View File

@ -39,6 +39,7 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
integer :: p integer :: p
double precision :: Ec double precision :: Ec
double precision :: flow
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision,allocatable :: Re_SigC(:) double precision,allocatable :: Re_SigC(:)
double precision,allocatable :: Im_SigC(:) double precision,allocatable :: Im_SigC(:)
@ -65,11 +66,14 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),Re_eHF(nOrb),Im_eHF(nOrb)) Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),Re_eHF(nOrb),Im_eHF(nOrb))
Re_eHF(:) = real(eHF(:)) Re_eHF(:) = real(eHF(:))
Im_eHF(:) = aimag(eHF(:)) Im_eHF(:) = aimag(eHF(:))
flow = 100d0
! Frequency-dependent second-order contribution ! Frequency-dependent second-order contribution
if(regularize) then
call complex_cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z) call complex_cRGF2_reg_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
else
call complex_cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
end if
Re_eGFlin(:) = Re_eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) Re_eGFlin(:) = Re_eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGFlin(:) = Im_eHF(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) Im_eGFlin(:) = Im_eHF(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
@ -84,7 +88,7 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*) write(*,*)
call complex_cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_eGFlin,Im_eGFlin,& call complex_cRGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_eGFlin,Im_eGFlin,&
Re_eHF,Im_eHF,Re_eGF,Im_eGF,Re_Z,Im_Z) Re_eHF,Im_eHF,Re_eGF,Im_eGF,Re_Z,Im_Z)
end if end if

13
src/GF/complex_cRGF2_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine complex_cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,& subroutine complex_cRGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z) ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z)
! Compute the graphical solution of the complex GF2 QP equation ! Compute the graphical solution of the complex GF2 QP equation
@ -9,6 +9,8 @@ subroutine complex_cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
! Input variables ! Input variables
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: flow
logical, intent(in) :: reg
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
@ -60,10 +62,13 @@ subroutine complex_cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt) do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt)
nIt = nIt + 1 nIt = nIt + 1
if(reg) then
call complex_cRGF_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,Re_w,Im_w,Re_eOld,Im_eOld,ERI,& call complex_cRGF_reg_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,Re_w,Im_w,Re_eOld,Im_eOld,ERI,&
Re_SigC,Im_SigC,Re_dSigC,Im_dSigC) Re_SigC,Im_SigC,Re_dSigC,Im_dSigC)
else
call complex_cRGF_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,Re_w,Im_w,Re_eOld,Im_eOld,ERI,&
Re_SigC,Im_SigC,Re_dSigC,Im_dSigC)
end if
Re_f = Re_w - Re_eHF(p) - Re_SigC Re_f = Re_w - Re_eHF(p) - Re_SigC
Im_f = Im_w - Im_eHF(p) - Im_SigC Im_f = Im_w - Im_eHF(p) - Im_SigC
Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2) Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)

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

@ -0,0 +1,98 @@
subroutine complex_cRGF2_reg_self_energy_diag(flow,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) :: flow
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 :: s
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
s = flow
! 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) &
*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
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)
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)&
*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
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)
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

2
src/GF/complex_cRGF2_self_energy_diag.f90
View File

@ -90,4 +90,6 @@ subroutine complex_cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,Re_e,Im_e,ERI,Re_
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
deallocate(Re_DS,Im_DS) deallocate(Re_DS,Im_DS)
call vecout(nBas,Re_SigC)
call vecout(nBas,Im_SigC)
end subroutine end subroutine

88
src/GF/complex_cRGF_reg_SigC_dSigC.f90 Normal file
View File

@ -0,0 +1,88 @@
subroutine complex_cRGF_reg_SigC_dSigC(s,p,eta,nBas,nC,nO,nV,nR,Re_w,Im_w,Re_e,Im_e,ERI,Re_SigC,Im_SigC,Re_DS,Im_DS)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: eta
double precision,intent(in) :: s
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)
double precision,intent(in) :: Re_w
double precision,intent(in) :: Im_w
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps
double precision :: eta_tilde
complex*16 :: num
complex*16 :: z_dummy
! 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
! Compute GF2 self-energy
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = Re_w + Re_e(a) - Re_e(i) - Re_e(j)
eta_tilde = eta - Im_w + 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)&
*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC = Re_SigC + real(z_dummy)
Im_SigC = Im_SigC + aimag(z_dummy)
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 = Re_DS + real(z_dummy)
Im_DS = Im_DS + aimag(z_dummy)
end do
end do
end do
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = Re_w + Re_e(i) - Re_e(a) - Re_e(b)
eta_tilde = eta + Im_w - 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)&
*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2)))
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
Re_SigC = Re_SigC + real(z_dummy)
Im_SigC = Im_SigC + aimag(z_dummy)
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 = Re_DS + real(z_dummy)
Im_DS = Im_DS + aimag(z_dummy)
end do
end do
end do
end subroutine

8
src/GF/complex_evRGF2.f90
View File

@ -42,6 +42,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
integer :: nSCF integer :: nSCF
integer :: n_diis integer :: n_diis
double precision :: Ec double precision :: Ec
double precision :: flow
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: Conv double precision :: Conv
double precision :: rcond double precision :: rcond
@ -93,6 +94,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
rcond = 0d0 rcond = 0d0
Re_Z(:) = 0d0 Re_Z(:) = 0d0
Im_Z(:) = 0d0 Im_Z(:) = 0d0
flow = 100d0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main SCF loop ! Main SCF loop
@ -103,8 +105,8 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
! Frequency-dependent second-order contribution ! Frequency-dependent second-order contribution
if(regularize) then if(regularize) then
write(*,*) "Implement regularization !!"
!call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z) call complex_cRGF2_reg_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,Re_eGF,Im_eGF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
else else
@ -123,7 +125,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*) write(*,*)
call complex_cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,& call complex_cRGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,&
ERI,Re_eOld,Im_eOld,Re_eOld,Im_eOld,Re_eGF,Im_eGF,Re_Z,Im_Z) ERI,Re_eOld,Im_eOld,Re_eOld,Im_eOld,Re_eGF,Im_eGF,Re_Z,Im_Z)
eGF = cmplx(Re_eGF,Im_eGF,kind=8) eGF = cmplx(Re_eGF,Im_eGF,kind=8)
end if end if