PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 12:43:48 +01:00
Code Releases Activity

merge self-energy and Z in UGTpp

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-31 09:23:43 +02:00
parent e449be4d18
commit 5620965a87
6 changed files with 163 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/GT/UG0T0pp.f90
View File

@ -92,17 +92,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Memory allocation ! Memory allocation
allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), &
Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), &
rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), &
Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), &
Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), &
rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), &
Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), &
Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), &
rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
SigT(nBas,nspin),Z(nBas,nspin), & SigT(nBas,nspin),Z(nBas,nspin),eG0T0(nBas,nspin))
eG0T0(nBas,nspin))
!---------------------------------------------- !----------------------------------------------
! alpha-beta block ! alpha-beta block
@ -110,15 +109,12 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
ispin = 1 ispin = 1
iblock = 3 iblock = 3
! iblock = 1
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin)
call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPab,Om1ab(:)) call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPab,Om1ab(:))
call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHab,Om2ab(:)) call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHab,Om2ab(:))
@ -134,9 +130,6 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPaa,Om1aa(:)) call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPaa,Om1aa(:))
call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHaa,Om2aa(:)) call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHaa,Om2aa(:))
@ -152,9 +145,6 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPbb,Om1bb(:)) call print_excitation_energies('ppRPA@HF (N+2)',iblock,nPbb,Om1bb(:))
call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHbb,Om2bb(:)) call print_excitation_energies('ppRPA@HF (N-2)',iblock,nHbb,Om2bb(:))
@ -162,10 +152,6 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy
!---------------------------------------------- !----------------------------------------------
EcGM = 0d0
SigT(:,:) = 0d0
Z(:,:) = 0d0
!alpha-beta block !alpha-beta block
iblock = 3 iblock = 3
@ -187,11 +173,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
rho1bb,X2bb,Y2bb,rho2bb) rho1bb,X2bb,Y2bb,rho2bb)
call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,Om1bb,& call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,Om1bb,&
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z)
call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,&
Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - Z(:,:)) Z(:,:) = 1d0/(1d0 - Z(:,:))
@ -199,13 +181,14 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Solve the quasi-particle equation ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
if(linearize) then if(linearize) then
eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
else else
eG0T0(:,:) = eHF(:,:) + SigT(:,:) write(*,*) 'Root search not yet implemented for UG0T0pp! Sorry.'
stop
end if end if

108
src/GT/UGTpp_renormalization_factor.f90
View File

@ -1,108 +0,0 @@
subroutine UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab, &
Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
! Compute renormalization factor of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC(nspin),nO(nspin),nV(nspin),nR(nspin)
integer,intent(in) :: nHaa,nHab,nHbb
integer,intent(in) :: nPaa,nPab,nPbb
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Om1aa(nPaa),Om1ab(nPab),Om1bb(nPbb)
double precision,intent(in) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab)
double precision,intent(in) :: rho1bb(nBas,nBas,nPbb)
double precision,intent(in) :: Om2aa(nHaa),Om2ab(nHab),Om2bb(nHbb)
double precision,intent(in) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab)
double precision,intent(in) :: rho2bb(nBas,nBas,nHbb)
! Local variables
integer :: i,a,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas,nspin)
!spin up part
! Occupied part of the T-matrix self-energy
do p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1)
do cd=1,nPaa
eps = e(p,1) + e(i,1) - Om1aa(cd)
Z(p,1) = Z(p,1) - rho1aa(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
do i=nC(1)+1,nO(1)
do cd=1,nPab
eps = e(p,1) + e(i,1) - Om1ab(cd)
Z(p,1) = Z(p,1) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
end do
end do
enddo
! Virtual part of the T-matrix self-energy
do p=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa
eps = e(p,1) + e(a,1) - Om2aa(kl)
Z(p,1) = Z(p,1) - rho2aa(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab
eps = e(p,1) + e(a,1) - Om2ab(kl)
Z(p,1) = Z(p,1) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
!spin down part
! Occupied part of the T-matrix self-energy
do p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2)
do cd=1,nPbb
eps = e(p,2) + e(i,2) - Om1bb(cd)
Z(p,2) = Z(p,2) - rho1bb(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
do i=nC(1)+1,nO(1)
do cd=1,nPab
eps = e(p,2) + e(i,2) - Om1ab(cd)
Z(p,2) = Z(p,2) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb
eps = e(p,2) + e(a,2) - Om2bb(kl)
Z(p,2) = Z(p,2) - rho2bb(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab
eps = e(p,2) + e(a,2) - Om2ab(kl)
Z(p,2) = Z(p,2) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo
enddo
enddo
end subroutine

111
src/GT/UGTpp_self_energy.f90
View File

@ -1,5 +1,5 @@
subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,& subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,&
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z)
! Compute the correlation part of the T-matrix self-energy ! Compute the correlation part of the T-matrix self-energy
@ -27,12 +27,19 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
! Local variables ! Local variables
integer :: i,j,a,b,p,q,cd,kl integer :: i,j,a,b,p,q,cd,kl
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(inout) :: EcGM(nspin) double precision,intent(inout) :: EcGM(nspin)
double precision,intent(inout) :: SigT(nBas,nBas,nspin) double precision,intent(inout) :: SigT(nBas,nBas,nspin)
double precision,intent(inout) :: Z(nBas,nspin)
! Initialization
EcGM(:) = 0d0
SigT(:,:,:) = 0d0
Z(:,:) = 0d0
!---------------------------------------------- !----------------------------------------------
! Occupied part of the T-matrix self-energy ! Occupied part of the T-matrix self-energy
@ -45,18 +52,22 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(p,1) + e(i,1) - Om1aa(cd) eps = e(p,1) + e(i,1) - Om1aa(cd)
SigT(p,q,1) = SigT(p,q,1) + rho1aa(p,i,cd)*rho1aa(q,i,cd)*eps/(eps**2 + eta**2) num = rho1aa(p,i,cd)*rho1aa(q,i,cd)
enddo SigT(p,q,1) = SigT(p,q,1) + num*eps/(eps**2 + eta**2)
enddo if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,1) + e(i,1) - Om1ab(cd) eps = e(p,1) + e(i,1) - Om1ab(cd)
SigT(p,q,1) = SigT(p,q,1) + rho1ab(p,i,cd)*rho1ab(q,i,cd)*eps/(eps**2 + eta**2) num = rho1ab(p,i,cd)*rho1ab(q,i,cd)
SigT(p,q,1) = SigT(p,q,1) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
enddo end do
!spin down part !spin down part
@ -65,18 +76,22 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPbb do cd=1,nPbb
eps = e(p,2) + e(i,2) - Om1bb(cd) eps = e(p,2) + e(i,2) - Om1bb(cd)
SigT(p,q,2) = SigT(p,q,2) + rho1bb(p,i,cd)*rho1bb(q,i,cd)*eps/(eps**2 + eta**2) num = rho1bb(p,i,cd)*rho1bb(q,i,cd)
enddo SigT(p,q,2) = SigT(p,q,2) + num*eps/(eps**2 + eta**2)
enddo if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,2) + e(i,2) - Om1ab(cd) eps = e(p,2) + e(i,2) - Om1ab(cd)
SigT(p,q,2) = SigT(p,q,2) + rho1ab(p,i,cd)*rho1ab(q,i,cd)*eps/(eps**2 + eta**2) num = rho1ab(p,i,cd)*rho1ab(q,i,cd)
SigT(p,q,2) = SigT(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
end do end do
end do end do
enddo end do
enddo end do
!---------------------------------------------- !----------------------------------------------
! Virtual part of the T-matrix self-energy ! Virtual part of the T-matrix self-energy
@ -89,18 +104,22 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa do kl=1,nHaa
eps = e(p,1) + e(a,1) - Om2aa(kl) eps = e(p,1) + e(a,1) - Om2aa(kl)
SigT(p,q,1) = SigT(p,q,1) + rho2aa(p,a,kl)*rho2aa(q,a,kl)*eps/(eps**2 + eta**2) num = rho2aa(p,a,kl)*rho2aa(q,a,kl)
enddo SigT(p,q,1) = SigT(p,q,1) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(p,1) + e(a,1) - Om2ab(kl) eps = e(p,1) + e(a,1) - Om2ab(kl)
SigT(p,q,1) = SigT(p,q,1) + rho2ab(p,a,kl)*rho2ab(q,a,kl)*eps/(eps**2 + eta**2) num = rho2ab(p,a,kl)*rho2ab(q,a,kl)
SigT(p,q,1) = SigT(p,q,1) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
enddo end do
!spin down part !spin down part
@ -109,18 +128,24 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(p,2) + e(a,2) - Om2bb(kl) eps = e(p,2) + e(a,2) - Om2bb(kl)
SigT(p,q,2) = SigT(p,q,2) + rho2bb(p,a,kl)*rho2bb(q,a,kl)*eps/(eps**2 + eta**2) num = rho2bb(p,a,kl)*rho2bb(q,a,kl)
enddo SigT(p,q,2) = SigT(p,q,2) + num*eps/(eps**2 + eta**2)
enddo if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab do kl=1,nHab
eps = e(p,2) + e(a,2) - Om2ab(kl) eps = e(p,2) + e(a,2) - Om2ab(kl)
SigT(p,q,2) = SigT(p,q,2) + rho2ab(p,a,kl)*rho2ab(q,a,kl)*eps/(eps**2 + eta**2) num = rho2ab(p,a,kl)*rho2ab(q,a,kl)
SigT(p,q,2) = SigT(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
end do end do
end do end do
enddo end do
enddo end do
Z(:,:) = 1d0/(1d0 - Z(:,:))
!---------------------------------------------- !----------------------------------------------
! Galitskii-Migdal correlation energy ! Galitskii-Migdal correlation energy
@ -133,9 +158,9 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do cd=1,nPaa do cd=1,nPaa
eps = e(i,1) + e(j,1) - Om1aa(cd) eps = e(i,1) + e(j,1) - Om1aa(cd)
EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
@ -151,18 +176,18 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do kl=1,nHaa do kl=1,nHaa
eps = e(a,1) + e(b,1) - Om2aa(kl) eps = e(a,1) + e(b,1) - Om2aa(kl)
EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(a,1) + e(b,1) - Om2ab(kl) eps = e(a,1) + e(b,1) - Om2ab(kl)
EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
! spin down part ! spin down part
@ -171,9 +196,9 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do cd=1,nPbb do cd=1,nPbb
eps = e(i,2) + e(j,2) - Om1bb(cd) eps = e(i,2) + e(j,2) - Om1bb(cd)
EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
@ -189,17 +214,17 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,
do kl=1,nHab do kl=1,nHab
eps = e(a,2) + e(b,2) - Om2ab(kl) eps = e(a,2) + e(b,2) - Om2ab(kl)
EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(a,2) + e(b,2) - Om2bb(kl) eps = e(a,2) + e(b,2) - Om2bb(kl)
EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
end subroutine end subroutine

111
src/GT/UGTpp_self_energy_diag.f90
View File

@ -1,5 +1,5 @@
subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,& subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,&
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z)
! Compute diagonal of the correlation part of the T-matrix self-energy ! Compute diagonal of the correlation part of the T-matrix self-energy
@ -27,52 +27,67 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,
! Local variables ! Local variables
integer :: i,j,a,b,p,cd,kl integer :: i,j,a,b,p,cd,kl
double precision :: eps double precision :: num,eps
! Output variables ! Output variables
double precision,intent(inout) :: EcGM(nspin) double precision,intent(inout) :: EcGM(nspin)
double precision,intent(inout) :: SigT(nBas,nspin) double precision,intent(inout) :: SigT(nBas,nspin)
double precision,intent(inout) :: Z(nBas,nspin)
! Initialization
EcGM(:) = 0d0
SigT(:,:) = 0d0
Z(:,:) = 0d0
!---------------------------------------------- !----------------------------------------------
! Occupied part of the T-matrix self-energy ! Occupied part of the T-matrix self-energy
!---------------------------------------------- !----------------------------------------------
!spin up part ! spin up part
do p=nC(1)+1,nBas-nR(1) do p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(p,1) + e(i,1) - Om1aa(cd) eps = e(p,1) + e(i,1) - Om1aa(cd)
SigT(p,1) = SigT(p,1) + rho1aa(p,i,cd)**2*eps/(eps**2 + eta**2) num = rho1aa(p,i,cd)**2
enddo SigT(p,1) = SigT(p,1) + num*eps/(eps**2 + eta**2)
enddo Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,1) + e(i,1) - Om1ab(cd) eps = e(p,1) + e(i,1) - Om1ab(cd)
SigT(p,1) = SigT(p,1) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) num = rho1ab(p,i,cd)**2
SigT(p,1) = SigT(p,1) + num*eps/(eps**2 + eta**2)
Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
!spin down part ! spin down part
do p=nC(2)+1,nBas-nR(2) do p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPbb do cd=1,nPbb
eps = e(p,2) + e(i,2) - Om1bb(cd) eps = e(p,2) + e(i,2) - Om1bb(cd)
SigT(p,2) = SigT(p,2) + rho1bb(p,i,cd)**2*eps/(eps**2 + eta**2) num = rho1bb(p,i,cd)**2
enddo SigT(p,2) = SigT(p,2) + num*eps/(eps**2 + eta**2)
enddo Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,2) + e(i,2) - Om1ab(cd) eps = e(p,2) + e(i,2) - Om1ab(cd)
SigT(p,2) = SigT(p,2) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) num = rho1ab(p,i,cd)**2
SigT(p,2) = SigT(p,2) + num*eps/(eps**2 + eta**2)
Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
!---------------------------------------------- !----------------------------------------------
! Virtual part of the T-matrix self-energy ! Virtual part of the T-matrix self-energy
@ -84,52 +99,60 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa do kl=1,nHaa
eps = e(p,1) + e(a,1) - Om2aa(kl) eps = e(p,1) + e(a,1) - Om2aa(kl)
SigT(p,1) = SigT(p,1) + rho2aa(p,a,kl)**2*eps/(eps**2 + eta**2) num = rho2aa(p,a,kl)**2
enddo SigT(p,1) = SigT(p,1) + num*eps/(eps**2 + eta**2)
Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(p,1) + e(a,1) - Om2ab(kl) eps = e(p,1) + e(a,1) - Om2ab(kl)
SigT(p,1) = SigT(p,1) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) num = rho2ab(p,a,kl)**2
SigT(p,1) = SigT(p,1) + num*eps/(eps**2 + eta**2)
Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
!spin down part ! spin down part
do p=nC(2)+1,nBas-nR(2) do p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(p,2) + e(a,2) - Om2bb(kl) eps = e(p,2) + e(a,2) - Om2bb(kl)
SigT(p,2) = SigT(p,2) + rho2bb(p,a,kl)**2*eps/(eps**2 + eta**2) num = rho2bb(p,a,kl)**2
enddo SigT(p,2) = SigT(p,2) + num*eps/(eps**2 + eta**2)
enddo Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab do kl=1,nHab
eps = e(p,2) + e(a,2) - Om2ab(kl) eps = e(p,2) + e(a,2) - Om2ab(kl)
SigT(p,2) = SigT(p,2) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) num = rho2ab(p,a,kl)**2
SigT(p,2) = SigT(p,2) + num*eps/(eps**2 + eta**2)
Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do end do
end do end do
enddo end do
Z(:,:) = 1d0/(1d0 - Z(:,:))
!---------------------------------------------- !----------------------------------------------
! Galitskii-Migdal correlation energy ! Galitskii-Migdal correlation energy
!---------------------------------------------- !----------------------------------------------
!spin up part ! spin up part
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do j=nC(1)+1,nO(1) do j=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(i,1) + e(j,1) - Om1aa(cd) eps = e(i,1) + e(j,1) - Om1aa(cd)
EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
@ -145,18 +168,18 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,
do kl=1,nHaa do kl=1,nHaa
eps = e(a,1) + e(b,1) - Om2aa(kl) eps = e(a,1) + e(b,1) - Om2aa(kl)
EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(a,1) + e(b,1) - Om2ab(kl) eps = e(a,1) + e(b,1) - Om2ab(kl)
EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
! spin down part ! spin down part
@ -165,9 +188,9 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,
do cd=1,nPbb do cd=1,nPbb
eps = e(i,2) + e(j,2) - Om1bb(cd) eps = e(i,2) + e(j,2) - Om1bb(cd)
EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
@ -183,17 +206,17 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,
do kl=1,nHab do kl=1,nHab
eps = e(a,2) + e(b,2) - Om2ab(kl) eps = e(a,2) + e(b,2) - Om2ab(kl)
EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(a,2) + e(b,2) - Om2bb(kl) eps = e(a,2) + e(b,2) - Om2bb(kl)
EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2)
enddo end do
enddo end do
enddo end do
end subroutine end subroutine

35
src/GT/evUGTpp.f90
View File

@ -98,17 +98,17 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Memory allocation ! Memory allocation
allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), &
Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), &
rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), &
Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), &
Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), &
rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), &
Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), &
Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), &
rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
SigT(nBas,nspin),Z(nBas,nspin), & SigT(nBas,nspin),Z(nBas,nspin),eGT(nBas,nspin), &
eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), & eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), &
e_diis(nBas,max_diis,nspin)) e_diis(nBas,max_diis,nspin))
!---------------------------------------------- !----------------------------------------------
@ -139,7 +139,6 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
ispin = 1 ispin = 1
iblock = 3 iblock = 3
! iblock = 1
! Compute linear response ! Compute linear response
@ -174,10 +173,6 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy
!---------------------------------------------- !----------------------------------------------
EcGM = 0d0
SigT(:,:) = 0d0
Z(:,:) = 0d0
!alpha-beta block !alpha-beta block
iblock = 3 iblock = 3
@ -199,13 +194,7 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
rho1bb,X2bb,Y2bb,rho2bb) rho1bb,X2bb,Y2bb,rho2bb)
call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,&
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z)
call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,&
Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - Z(:,:))
!---------------------------------------------- !----------------------------------------------
! Solve the quasi-particle equation ! Solve the quasi-particle equation

14
src/GT/qsUGTpp.f90
View File

@ -248,10 +248,6 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy
!---------------------------------------------- !----------------------------------------------
EcGM = 0d0
SigT(:,:,:) = 0d0
Z(:,:) = 0d0
!alpha-beta block !alpha-beta block
iblock = 3 iblock = 3
@ -273,13 +269,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
rho1bb,X2bb,Y2bb,rho2bb) rho1bb,X2bb,Y2bb,rho2bb)
call UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& call UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,&
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z)
call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,&
Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - Z(:,:))
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
@ -296,7 +286,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
do ispin=1,nspin do ispin=1,nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) & F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) &
+ SigTp(:,:,ispin) + SigTp(:,:,ispin)
end do end do
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria