From 541181594d4ae4b9f3ddca07091a60a68a6d0bf7 Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Wed, 26 Mar 2025 12:32:36 +0100 Subject: [PATCH] root search for cG0F2 --- src/GF/cRG0F2.f90 | 2 +- src/GF/cRGF2_Im_SigC.f90 | 59 ++++++++++++++++++++++++ src/GF/cRGF2_Im_dSigC.f90 | 56 ++++++++++++++++++++++ src/GF/cRGF2_QP_graph.f90 | 97 +++++++++++++++++++++++++++++++++++++++ src/GF/cRGF2_Re_SigC.f90 | 58 +++++++++++++++++++++++ src/GF/cRGF2_Re_dSigC.f90 | 59 ++++++++++++++++++++++++ 6 files changed, 330 insertions(+), 1 deletion(-) create mode 100644 src/GF/cRGF2_Im_SigC.f90 create mode 100644 src/GF/cRGF2_Im_dSigC.f90 create mode 100644 src/GF/cRGF2_QP_graph.f90 create mode 100644 src/GF/cRGF2_Re_SigC.f90 create mode 100644 src/GF/cRGF2_Re_dSigC.f90 diff --git a/src/GF/cRG0F2.f90 b/src/GF/cRG0F2.f90 index 4e05d89..e9e1f51 100644 --- a/src/GF/cRG0F2.f90 +++ b/src/GF/cRG0F2.f90 @@ -93,7 +93,7 @@ subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize write(*,*) ' *** Quasiparticle energies obtained by root search *** NOT IMPLEMEMTED YET' write(*,*) - !call RGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z) + 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 diff --git a/src/GF/cRGF2_Im_SigC.f90 b/src/GF/cRGF2_Im_SigC.f90 new file mode 100644 index 0000000..8e9338e --- /dev/null +++ b/src/GF/cRGF2_Im_SigC.f90 @@ -0,0 +1,59 @@ +double precision function cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: Re_w + double precision,intent(in) :: Im_w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + double precision :: eta_tilde + double precision :: num + +! Initialize + + cRGF2_Im_SigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = Re_w + eHF(a) - eHF(i) - eHF(j) + eta_tilde = eta - Im_w + e_cap(i) + e_cap(a) - e_cap(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + cRGF2_Im_SigC = cRGF2_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2) + + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = Re_w + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) + cRGF2_Im_SigC = cRGF2_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2) + end do + end do + end do + +end function diff --git a/src/GF/cRGF2_Im_dSigC.f90 b/src/GF/cRGF2_Im_dSigC.f90 new file mode 100644 index 0000000..1678817 --- /dev/null +++ b/src/GF/cRGF2_Im_dSigC.f90 @@ -0,0 +1,56 @@ +double precision function cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: Re_w + double precision,intent(in) :: Im_w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + double precision :: eta_tilde + double precision :: num + + ! Initialize + + cRGF2_Im_dSigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + eps = Re_w + eHF(a) - eHF(i) - eHF(j) + eta_tilde = eta - Im_w + e_cap(i) + e_cap(a) - e_cap(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + cRGF2_Im_dSigC = cRGF2_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + eps = Re_w + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) + cRGF2_Im_dSigC = cRGF2_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 + end do + end do + end do + +end function diff --git a/src/GF/cRGF2_QP_graph.f90 b/src/GF/cRGF2_QP_graph.f90 new file mode 100644 index 0000000..2f2a6ff --- /dev/null +++ b/src/GF/cRGF2_QP_graph.f90 @@ -0,0 +1,97 @@ +subroutine cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z) + +! Compute the graphical solution of the GF2 QP equation + + 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) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: Re_eGFlin(nBas) + double precision,intent(in) :: Im_eGFlin(nBas) + double precision,intent(in) :: Re_eOld(nBas) + double precision,intent(in) :: Im_eOld(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision,external :: cRGF2_Re_SigC,cRGF2_Im_SigC,cRGF2_Re_dSigC,cRGF2_Im_dSigC + double precision :: Re_SigC,Im_SigC,Re_dSigC,Im_dSigC + double precision :: Re_f,Im_f,Re_df,Im_df + double precision :: Re_w,Im_w + +! Output variables + + double precision,intent(out) :: Re_eGF(nBas),Im_eGF(nBas) + double precision,intent(out) :: Re_Z(nBas),Im_Z(nBas) + +! Run Newton's algorithm to find the root + + write(*,*)'-----------------------------------------------------' + write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','Re(e_GFlin) (eV)','Re(e_GF) (eV)','Re(Z)' + write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','Im(e_GFlin) (eV)','Im(e_GF) (eV)','Im(Z)' + write(*,*)'-----------------------------------------------------' + + do p=nC+1,nBas-nR + + Re_w = Re_eGFlin(p) + Im_w = Im_eGFlin(p) + nIt = 0 + Re_f = 1d0 + Im_f = 0d0 + + do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + + Re_SigC = cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) + Im_SigC = cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) + Re_dSigC = cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) + Im_dSigC = cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) + + Re_f = Re_w - eHF(p) - Re_SigC + Im_f = Im_w - e_cap(p) - Im_SigC + Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2) + Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2) + + Re_w = Re_w - Re_df*Re_f + Im_df*Im_f + Im_w = Im_w - Re_f*Im_df - Re_df*Im_f + + end do + + if(nIt == maxIt) then + + Re_eGF(p) = Re_eGFlin(p) + Im_eGF(p) = Im_eGFlin(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p),'Cvg Failed!' + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p),'Cvg Failed!' + + else + + Re_eGF(p) = Re_w + Im_eGF(p) = Im_w + Re_Z(p) = Re_df + Im_Z(p) = Im_df + + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p) + + write(*,*)'-----------------------------------------------------' + end if + + end do + +end subroutine diff --git a/src/GF/cRGF2_Re_SigC.f90 b/src/GF/cRGF2_Re_SigC.f90 new file mode 100644 index 0000000..11329f1 --- /dev/null +++ b/src/GF/cRGF2_Re_SigC.f90 @@ -0,0 +1,58 @@ +double precision function cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: Re_w + double precision,intent(in) :: Im_w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + double precision :: eta_tilde + double precision :: num + +! Initialize + + cRGF2_Re_SigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = Re_w + eHF(a) - eHF(i) - eHF(j) + eta_tilde = eta - Im_w + e_cap(i) + e_cap(a) - e_cap(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2) + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = Re_w + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) + cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2) + end do + end do + end do + +end function diff --git a/src/GF/cRGF2_Re_dSigC.f90 b/src/GF/cRGF2_Re_dSigC.f90 new file mode 100644 index 0000000..2d74c8f --- /dev/null +++ b/src/GF/cRGF2_Re_dSigC.f90 @@ -0,0 +1,59 @@ +double precision function cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: Re_w + double precision,intent(in) :: Im_w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + double precision :: eta_tilde + double precision :: num +! Initialize + + cRGF2_Re_dSigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = Re_w + eHF(a) - eHF(i) - eHF(j) + eta_tilde = eta - Im_w + e_cap(i) + e_cap(a) - e_cap(j) + cRGF2_Re_dSigC = cRGF2_Re_dSigC -& + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 + + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = Re_w + eHF(i) - eHF(a) - eHF(b) + eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) + cRGF2_Re_dSigC = cRGF2_Re_dSigC -& + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 + + end do + end do + end do + +end function