From 4334ce3ee6aa889de65dfd2f69942f1dcc83de6b Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 30 Jun 2023 19:17:35 +0200 Subject: [PATCH] missing files --- src/GF/QP_graph_GF2.f90 | 76 +++++++++++++++++++++++++++++++++++++++++ src/GF/SigmaC_GF2.f90 | 47 +++++++++++++++++++++++++ src/GF/dSigmaC_GF2.f90 | 47 +++++++++++++++++++++++++ src/GT/QP_graph_GT.f90 | 76 +++++++++++++++++++++++++++++++++++++++++ src/GT/SigmaC_GT.f90 | 55 +++++++++++++++++++++++++++++ src/GT/dSigmaC_GT.f90 | 54 +++++++++++++++++++++++++++++ 6 files changed, 355 insertions(+) create mode 100644 src/GF/QP_graph_GF2.f90 create mode 100644 src/GF/SigmaC_GF2.f90 create mode 100644 src/GF/dSigmaC_GF2.f90 create mode 100644 src/GT/QP_graph_GT.f90 create mode 100644 src/GT/SigmaC_GT.f90 create mode 100644 src/GT/dSigmaC_GT.f90 diff --git a/src/GF/QP_graph_GF2.f90 b/src/GF/QP_graph_GF2.f90 new file mode 100644 index 0000000..1017505 --- /dev/null +++ b/src/GF/QP_graph_GF2.f90 @@ -0,0 +1,76 @@ +subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2) + +! 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,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGF2lin(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 :: SigmaC_GF2,dSigmaC_GF2 + double precision :: sigC,dsigC + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGF2(nBas) + + +! Run Newton's algorithm to find the root + + do p=nC+1,nBas-nR + + write(*,*) '-----------------' + write(*,'(A10,I3)') 'Orbital ',p + write(*,*) '-----------------' + + w = eGF2lin(p) + nIt = 0 + f = 1d0 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + sigC = SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI) + dsigC = dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI) + f = w - eHF(p) - sigC + df = 1d0 - dsigC + + w = w - f/df + + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC + + + end do + + if(nIt == maxIt) then + + write(*,*) 'Newton root search has not converged!' + + else + + eGF2(p) = w + + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGF2(p)*HaToeV + write(*,*) + + end if + +end do + +end subroutine QP_graph_GF2 diff --git a/src/GF/SigmaC_GF2.f90 b/src/GF/SigmaC_GF2.f90 new file mode 100644 index 0000000..8f099fb --- /dev/null +++ b/src/GF/SigmaC_GF2.f90 @@ -0,0 +1,47 @@ +double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,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) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + + SigmaC_GF2 = 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 = w + eHF(a) - eHF(i) - eHF(j) + SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**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 = w + eHF(i) - eHF(a) - eHF(b) + SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2) + + end do + end do + end do + +end function SigmaC_GF2 diff --git a/src/GF/dSigmaC_GF2.f90 b/src/GF/dSigmaC_GF2.f90 new file mode 100644 index 0000000..4400584 --- /dev/null +++ b/src/GF/dSigmaC_GF2.f90 @@ -0,0 +1,47 @@ +double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,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) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + + + + ! Occupied part of the correlation self-energy + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = w + eHF(a) - eHF(i) - eHF(j) + dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**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 = w + eHF(i) - eHF(a) - eHF(b) + dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + +end function dSigmaC_GF2 diff --git a/src/GT/QP_graph_GT.f90 b/src/GT/QP_graph_GT.f90 new file mode 100644 index 0000000..75994bc --- /dev/null +++ b/src/GT/QP_graph_GT.f90 @@ -0,0 +1,76 @@ +subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,eGTlin,eGT) + +! Iput variables + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + + double precision,intent(in) :: eGTlin(nBas) + +! Local variables + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision,external :: SigmaC_GT,dSigmaC_GT + double precision :: sigC,dsigC + double precision :: f,df + double precision :: w + +! Output variables + double precision,intent(out) :: eGT(nBas) + + sigC = 0d0 + dsigC = 0d0 + +! Run Newton's algorithm to find the root + do p=nC+1,nBas-nR + + write(*,*) '-----------------' + write(*,'(A10,I3)') 'Orbital ',p + write(*,*) '-----------------' + + w = eGTlin(p) + write(*,*) 'HERE', eGTlin(p), eHF(p) + nIt = 0 + f = 1d0 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + sigC = SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2) + dsigC = dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2) + write (*,*) sigC + f = w - eHF(p) - sigC + df = 1d0 - dsigC + + w = w - f/df + + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC + + end do + + if(nIt == maxIt) then + write(*,*) 'Newton root search has not converged!' + else + eGT(p) = w + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV + write(*,*) + end if + + end do + +end subroutine QP_graph_GT diff --git a/src/GT/SigmaC_GT.f90 b/src/GT/SigmaC_GT.f90 new file mode 100644 index 0000000..770aee0 --- /dev/null +++ b/src/GT/SigmaC_GT.f90 @@ -0,0 +1,55 @@ +double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2) + +! 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) :: w + 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 + integer,intent(in) :: nOO,nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,a,cd,kl + double precision :: eps + +! Initialize + + SigmaC_GT = 0d0 + +!---------------------------------------------- +! Occupied part of the T-matrix self-energy +!---------------------------------------------- + do i=nC+1,nO + do cd=1,nVV + eps = w + e(i) - Omega1(cd) + SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + enddo + enddo + write (*,*) SigmaC_GT +!---------------------------------------------- +! Virtual part of the T-matrix self-energy +!---------------------------------------------- + do a=nO+1,nBas-nR + do kl=1,nOO + eps = w + e(a) - Omega2(kl) + SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + enddo + enddo + write (*,*) SigmaC_GT + +end function SigmaC_GT diff --git a/src/GT/dSigmaC_GT.f90 b/src/GT/dSigmaC_GT.f90 new file mode 100644 index 0000000..0e73b11 --- /dev/null +++ b/src/GT/dSigmaC_GT.f90 @@ -0,0 +1,54 @@ +double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2) + +! 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) :: w + 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 + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,a,cd,kl + double precision :: eps + +! Initialize + + dSigmaC_GT = 0d0 + +!---------------------------------------------- +! Occupied part of the T-matrix self-energy +!---------------------------------------------- + do i=nC+1,nO + do cd=1,nVV + eps = w + e(i) - Omega1(cd) + dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo +!---------------------------------------------- +! Virtual part of the T-matrix self-energy +!---------------------------------------------- + do a=nO+1,nBas-nR + do kl=1,nOO + eps = w + e(a) - Omega2(kl) + dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo + +end function dSigmaC_GT