From 02f7a03385f99ed7b2f3816b5e2d493783f200fb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 17 Dec 2021 13:36:26 +0100 Subject: [PATCH] regularized GF2 --- include/parameters.h | 1 - input/methods | 2 +- src/GF/G0F2.f90 | 10 +- src/GF/UG0F2.f90 | 10 +- src/GF/evGF2.f90 | 10 +- src/GF/evUGF2.f90 | 10 +- src/GF/qsGF2.f90 | 10 +- src/GF/qsUGF2.f90 | 10 +- src/GF/regularized_self_energy_GF2.f90 | 92 +++++++ src/GF/regularized_self_energy_GF2_diag.f90 | 88 +++++++ ...restricted_regularized_self_energy_GF2.f90 | 236 ++++++++++++++++++ ...icted_regularized_self_energy_GF2_diag.f90 | 231 +++++++++++++++++ 12 files changed, 702 insertions(+), 8 deletions(-) create mode 100644 src/GF/regularized_self_energy_GF2.f90 create mode 100644 src/GF/regularized_self_energy_GF2_diag.f90 create mode 100644 src/GF/unrestricted_regularized_self_energy_GF2.f90 create mode 100644 src/GF/unrestricted_regularized_self_energy_GF2_diag.f90 diff --git a/include/parameters.h b/include/parameters.h index 2644d6a..dbf9364 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -20,4 +20,3 @@ double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: CxLSDA = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - diff --git a/input/methods b/input/methods index df2ad59..196855c 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F T F F F + T F F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index 2e98066..391ed99 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -59,7 +59,15 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize ! Frequency-dependent second-order contribution - call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) + if(regularize) then + + call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) + + else + + call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) + + end if if(linearize) then diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 index be7d4c1..8c81da8 100644 --- a/src/GF/UG0F2.f90 +++ b/src/GF/UG0F2.f90 @@ -80,7 +80,15 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta, ! Compute self-energy ! !---------------------! - call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z) + + else + + call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z) + + end if !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index f55ae32..9466f53 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -79,7 +79,15 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, ! Frequency-dependent second-order contribution - call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + if(regularize) then + + call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + + else + + call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + + end if if(linearize) then diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 07749af..3d7c844 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -113,7 +113,15 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Compute self-energy and renormalization factor ! !------------------------------------------------! - call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + + else + + call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + + end if !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 3019091..6a5a395 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -145,7 +145,15 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ! Compute self-energy and renormalization factor - call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + if(regularize) then + + call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + + else + + call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + + end if ! Make correlation self-energy Hermitian and transform it back to AO basis diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 71d3c97..e94a344 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -179,7 +179,15 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Compute self-energy and renormalization factor ! !------------------------------------------------! - call unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + + else + + call unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + + end if ! Make correlation self-energy Hermitian and transform it back to AO basis diff --git a/src/GF/regularized_self_energy_GF2.f90 b/src/GF/regularized_self_energy_GF2.f90 new file mode 100644 index 0000000..e4b0475 --- /dev/null +++ b/src/GF/regularized_self_energy_GF2.f90 @@ -0,0 +1,92 @@ +subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + +! Compute GF2 self-energy and its renormalization factor + + 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) :: eGF2(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + integer :: p,q + double precision :: eps + double precision :: num + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: Z(nBas) + +! Initialize + + SigC(:,:) = 0d0 + Z(:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +!----------------------------------------------------! +! Compute GF2 self-energy and renormalization factor ! +!----------------------------------------------------! + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q) = SigC(p,q) + num*fk + if(p == q) Z(p) = Z(p) - num*dfk + + end do + end do + end do + end do + end do + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q) = SigC(p,q) + num*fk + if(p == q) Z(p) = Z(p) - num*dfk + + end do + end do + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine regularized_self_energy_GF2 diff --git a/src/GF/regularized_self_energy_GF2_diag.f90 b/src/GF/regularized_self_energy_GF2_diag.f90 new file mode 100644 index 0000000..80693e8 --- /dev/null +++ b/src/GF/regularized_self_energy_GF2_diag.f90 @@ -0,0 +1,88 @@ +subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,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,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGF2(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + integer :: p + double precision :: eps + double precision :: num + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: Z(nBas) + +! Initialize + + SigC(:) = 0d0 + Z(:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +!----------------------------------------------------! +! Compute GF2 self-energy and renormalization factor ! +!----------------------------------------------------! + + do p=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p) = SigC(p) + num*fk + Z(p) = Z(p) - num*dfk + + 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 = eGF2(p) + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p) = SigC(p) + num*fk + Z(p) = Z(p) - num*dfk + + end do + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine regularized_self_energy_GF2_diag diff --git a/src/GF/unrestricted_regularized_self_energy_GF2.f90 b/src/GF/unrestricted_regularized_self_energy_GF2.f90 new file mode 100644 index 0000000..dcae0e4 --- /dev/null +++ b/src/GF/unrestricted_regularized_self_energy_GF2.f90 @@ -0,0 +1,236 @@ +subroutine unrestricted_regularized_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,SigC,Z) + +! Perform unrestricted GF2 self-energy and its renormalization factor + + 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) + double precision,intent(in) :: eta + double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: eGF2(nBas,nspin) + +! Local variables + + integer :: p,q + integer :: i,j,a,b + double precision :: eps,num + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + +!---------------------! +! Compute self-energy | +!---------------------! + + SigC(:,:,:) = 0d0 + Z(:,:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + + !----------------! + ! Spin-up sector + !----------------! + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + + ! Addition part: aa + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + + eps = eGF2(p,1) + eHF(i,1) - eHF(a,1) - eHF(b,1) + num = ERI_aa(i,q,a,b)*ERI_aa(a,b,i,p) & + - ERI_aa(i,q,a,b)*ERI_aa(a,b,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,1) = SigC(p,q,1) + num*fk + if(p == q) Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Addition part: ab + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do b=nO(1)+1,nBas-nR(1) + + eps = eGF2(p,1) + eHF(i,2) - eHF(a,2) - eHF(b,1) + num = ERI_ab(q,i,b,a)*ERI_ab(b,a,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,1) = SigC(p,q,1) + num*fk + if(p == q) Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Removal part: aa + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do j=nC(1)+1,nO(1) + + eps = eGF2(p,1) + eHF(a,1) - eHF(i,1) - eHF(j,1) + num = ERI_aa(a,q,i,j)*ERI_aa(i,j,a,p) & + - ERI_aa(a,q,i,j)*ERI_aa(i,j,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,1) = SigC(p,q,1) + num*fk + if(p == q) Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Removal part: ab + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do j=nC(1)+1,nO(1) + + eps = eGF2(p,1) + eHF(a,2) - eHF(i,2) - eHF(j,1) + num = ERI_ab(q,a,j,i)*ERI_ab(j,i,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,1) = SigC(p,q,1) + num*fk + if(p == q) Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + enddo + enddo + + !------------------! + ! Spin-down sector ! + !------------------! + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + + ! Addition part: bb + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do b=nO(2)+1,nBas-nR(2) + + eps = eGF2(p,2) + eHF(i,2) - eHF(a,2) - eHF(b,2) + num = ERI_bb(i,q,a,b)*ERI_bb(a,b,i,p) & + - ERI_bb(i,q,a,b)*ERI_bb(a,b,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + enddo + enddo + + ! Addition part: ab + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do b=nO(2)+1,nBas-nR(2) + + eps = eGF2(p,2) + eHF(i,1) - eHF(a,1) - eHF(b,2) + num = ERI_ab(i,q,a,b)*ERI_ab(a,b,i,p) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,2) = SigC(p,q,2) + num*fk + if(p == q) Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + ! Removal part: bb + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do j=nC(2)+1,nO(2) + + eps = eGF2(p,2) + eHF(a,2) - eHF(i,2) - eHF(j,2) + num = ERI_bb(a,q,i,j)*ERI_bb(i,j,a,p) & + - ERI_bb(a,q,i,j)*ERI_bb(i,j,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,2) = SigC(p,q,2) + num*fk + if(p == q) Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + ! Removal part: ab + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do j=nC(2)+1,nO(2) + + eps = eGF2(p,2) + eHF(a,1) - eHF(i,1) - eHF(j,2) + num = ERI_ab(a,q,i,j)*ERI_ab(i,j,a,p) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,q,2) = SigC(p,q,2) + num*fk + if(p == q) Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + enddo + enddo + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +end subroutine unrestricted_regularized_self_energy_GF2 diff --git a/src/GF/unrestricted_regularized_self_energy_GF2_diag.f90 b/src/GF/unrestricted_regularized_self_energy_GF2_diag.f90 new file mode 100644 index 0000000..6818d0e --- /dev/null +++ b/src/GF/unrestricted_regularized_self_energy_GF2_diag.f90 @@ -0,0 +1,231 @@ +subroutine unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,SigC,Z) + +! Perform unrestricted GF2 self-energy and its renormalization factor + + 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) + double precision,intent(in) :: eta + double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: eGF2(nBas,nspin) + +! Local variables + + integer :: p + integer :: i,j,a,b + double precision :: eps,num + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + +!---------------------! +! Compute self-energy | +!---------------------! + + SigC(:,:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + + !----------------! + ! Spin-up sector + !----------------! + + do p=nC(1)+1,nBas-nR(1) + + ! Addition part: aa + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + + eps = eGF2(p,1) + eHF(i,1) - eHF(a,1) - eHF(b,1) + num = ERI_aa(i,p,a,b)*ERI_aa(a,b,i,p) & + - ERI_aa(i,p,a,b)*ERI_aa(a,b,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,1) = SigC(p,1) + num*fk + Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Addition part: ab + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do b=nO(1)+1,nBas-nR(1) + + eps = eGF2(p,1) + eHF(i,2) - eHF(a,2) - eHF(b,1) + num = ERI_ab(p,i,b,a)*ERI_ab(b,a,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,1) = SigC(p,1) + num*fk + Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Removal part: aa + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do j=nC(1)+1,nO(1) + + eps = eGF2(p,1) + eHF(a,1) - eHF(i,1) - eHF(j,1) + num = ERI_aa(a,p,i,j)*ERI_aa(i,j,a,p) & + - ERI_aa(a,p,i,j)*ERI_aa(i,j,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,1) = SigC(p,1) + num*fk + Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + ! Removal part: ab + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do j=nC(1)+1,nO(1) + + eps = eGF2(p,1) + eHF(a,2) - eHF(i,2) - eHF(j,1) + num = ERI_ab(p,a,j,i)*ERI_ab(j,i,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,1) = SigC(p,1) + num*fk + Z(p,1) = Z(p,1) - num*dfk + + enddo + enddo + enddo + + enddo + + !------------------! + ! Spin-down sector ! + !------------------! + + do p=nC(2)+1,nBas-nR(2) + + ! Addition part: bb + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do b=nO(2)+1,nBas-nR(2) + + eps = eGF2(p,2) + eHF(i,2) - eHF(a,2) - eHF(b,2) + num = ERI_bb(i,p,a,b)*ERI_bb(a,b,i,p) & + - ERI_bb(i,p,a,b)*ERI_bb(a,b,p,i) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,2) = SigC(p,2) + num*fk + Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + ! Addition part: ab + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do b=nO(2)+1,nBas-nR(2) + + eps = eGF2(p,2) + eHF(i,1) - eHF(a,1) - eHF(b,2) + num = ERI_ab(i,p,a,b)*ERI_ab(a,b,i,p) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,2) = SigC(p,2) + num*fk + Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + ! Removal part: bb + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do j=nC(2)+1,nO(2) + + eps = eGF2(p,2) + eHF(a,2) - eHF(i,2) - eHF(j,2) + num = ERI_bb(a,p,i,j)*ERI_bb(i,j,a,p) & + - ERI_bb(a,p,i,j)*ERI_bb(i,j,p,a) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,2) = SigC(p,2) + num*fk + Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + ! Removal part: ab + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do j=nC(2)+1,nO(2) + + eps = eGF2(p,2) + eHF(a,1) - eHF(i,1) - eHF(j,2) + num = ERI_ab(a,p,i,j)*ERI_ab(i,j,a,p) + + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + + SigC(p,2) = SigC(p,2) + num*fk + Z(p,2) = Z(p,2) - num*dfk + + enddo + enddo + enddo + + enddo + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +end subroutine unrestricted_regularized_self_energy_GF2_diag