From 5620965a87b9ea9afee7c2f38f7da984b2372cbb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 31 Jul 2023 09:23:43 +0200 Subject: [PATCH] merge self-energy and Z in UGTpp --- src/GT/UG0T0pp.f90 | 45 +++------- src/GT/UGTpp_renormalization_factor.f90 | 108 ----------------------- src/GT/UGTpp_self_energy.f90 | 111 +++++++++++++++--------- src/GT/UGTpp_self_energy_diag.f90 | 111 ++++++++++++++---------- src/GT/evUGTpp.f90 | 35 +++----- src/GT/qsUGTpp.f90 | 14 +-- 6 files changed, 163 insertions(+), 261 deletions(-) delete mode 100644 src/GT/UGTpp_renormalization_factor.f90 diff --git a/src/GT/UG0T0pp.f90 b/src/GT/UG0T0pp.f90 index 09c63cc..a34d926 100644 --- a/src/GT/UG0T0pp.f90 +++ b/src/GT/UG0T0pp.f90 @@ -92,17 +92,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! Memory allocation - allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & - Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & - Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & - Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & - Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & - Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & - SigT(nBas,nspin),Z(nBas,nspin), & - eG0T0(nBas,nspin)) + SigT(nBas,nspin),Z(nBas,nspin),eG0T0(nBas,nspin)) !---------------------------------------------- ! alpha-beta block @@ -110,15 +109,12 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ispin = 1 iblock = 3 -! iblock = 1 ! Compute linear response 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)) -! 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,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, & 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,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, & 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,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 !---------------------------------------------- - EcGM = 0d0 - SigT(:,:) = 0d0 - Z(:,:) = 0d0 - !alpha-beta block iblock = 3 @@ -187,11 +173,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & 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,& - rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) - - 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) - + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,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 !---------------------------------------------- - if(linearize) then + if(linearize) then + + eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:) - eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:) - else - - eG0T0(:,:) = eHF(:,:) + SigT(:,:) + + write(*,*) 'Root search not yet implemented for UG0T0pp! Sorry.' + stop end if diff --git a/src/GT/UGTpp_renormalization_factor.f90 b/src/GT/UGTpp_renormalization_factor.f90 deleted file mode 100644 index a66e1fc..0000000 --- a/src/GT/UGTpp_renormalization_factor.f90 +++ /dev/null @@ -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 diff --git a/src/GT/UGTpp_self_energy.f90 b/src/GT/UGTpp_self_energy.f90 index ba2fd06..94da65b 100644 --- a/src/GT/UGTpp_self_energy.f90 +++ b/src/GT/UGTpp_self_energy.f90 @@ -1,5 +1,5 @@ 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 @@ -27,12 +27,19 @@ subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb, ! Local variables integer :: i,j,a,b,p,q,cd,kl - double precision :: eps + double precision :: num,eps ! Output variables double precision,intent(inout) :: EcGM(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 @@ -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 cd=1,nPaa 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) - enddo - enddo + num = rho1aa(p,i,cd)*rho1aa(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 do i=nC(2)+1,nO(2) do cd=1,nPab 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 - enddo - enddo + end do + end do !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 cd=1,nPbb 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) - enddo - enddo + num = rho1bb(p,i,cd)*rho1bb(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 do i=nC(2)+1,nO(2) do cd=1,nPab 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 - enddo - enddo + end do + end do !---------------------------------------------- ! 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 kl=1,nHaa 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) - enddo + num = rho2aa(p,a,kl)*rho2aa(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 do a=nO(1)+1,nBas-nR(1) do kl=1,nHab 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 - enddo - enddo + end do + end do !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 kl=1,nHbb 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) - enddo - enddo + num = rho2bb(p,a,kl)*rho2bb(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 do a=nO(2)+1,nBas-nR(2) do kl=1,nHab 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 - enddo - enddo + end do + end do + + Z(:,:) = 1d0/(1d0 - Z(:,:)) !---------------------------------------------- ! 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 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) - enddo - enddo - enddo + end do + end do + end do do i=nC(1)+1,nO(1) 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 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) - enddo - enddo - enddo + end do + end do + end do do a=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1) do kl=1,nHab 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) - enddo - enddo - enddo + end do + end do + end do ! 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 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) - enddo - enddo - enddo + end do + end do + end do do i=nC(1)+1,nO(1) 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 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) - enddo - enddo - enddo + end do + end do + end do do a=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2) do kl=1,nHbb 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) - enddo - enddo - enddo + end do + end do + end do end subroutine diff --git a/src/GT/UGTpp_self_energy_diag.f90 b/src/GT/UGTpp_self_energy_diag.f90 index 9582b13..46539bb 100644 --- a/src/GT/UGTpp_self_energy_diag.f90 +++ b/src/GT/UGTpp_self_energy_diag.f90 @@ -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,& - 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 @@ -27,52 +27,67 @@ subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab, ! Local variables integer :: i,j,a,b,p,cd,kl - double precision :: eps + double precision :: num,eps ! Output variables double precision,intent(inout) :: EcGM(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 !---------------------------------------------- -!spin up part +! spin up part 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) - SigT(p,1) = SigT(p,1) + rho1aa(p,i,cd)**2*eps/(eps**2 + eta**2) - enddo - enddo + num = rho1aa(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 do i=nC(2)+1,nO(2) do cd=1,nPab 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 - enddo + end do -!spin down part +! spin down part 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) - SigT(p,2) = SigT(p,2) + rho1bb(p,i,cd)**2*eps/(eps**2 + eta**2) - enddo - enddo + num = rho1bb(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 do i=nC(2)+1,nO(2) do cd=1,nPab 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 - enddo + end do !---------------------------------------------- ! 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 kl=1,nHaa 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) - enddo + num = rho2aa(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 do a=nO(1)+1,nBas-nR(1) do kl=1,nHab 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 - enddo + end do -!spin down part +! spin down part 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) - SigT(p,2) = SigT(p,2) + rho2bb(p,a,kl)**2*eps/(eps**2 + eta**2) - enddo - enddo + num = rho2bb(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 do a=nO(2)+1,nBas-nR(2) do kl=1,nHab 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 - enddo - + end do + Z(:,:) = 1d0/(1d0 - Z(:,:)) !---------------------------------------------- ! Galitskii-Migdal correlation energy !---------------------------------------------- -!spin up part +! spin up part do i=nC(1)+1,nO(1) do j=nC(1)+1,nO(1) do cd=1,nPaa 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) - enddo - enddo - enddo + end do + end do + end do do i=nC(1)+1,nO(1) 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 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) - enddo - enddo - enddo + end do + end do + end do do a=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1) do kl=1,nHab 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) - enddo - enddo - enddo + end do + end do + end do ! 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 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) - enddo - enddo - enddo + end do + end do + end do do i=nC(1)+1,nO(1) 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 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) - enddo - enddo - enddo + end do + end do + end do do a=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2) do kl=1,nHbb 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) - enddo - enddo - enddo + end do + end do + end do end subroutine diff --git a/src/GT/evUGTpp.f90 b/src/GT/evUGTpp.f90 index 23b7e17..5217ace 100644 --- a/src/GT/evUGTpp.f90 +++ b/src/GT/evUGTpp.f90 @@ -98,17 +98,17 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Memory allocation - allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & - Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & - rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & - Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & - Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & - rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & - Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & - Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & - rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & - SigT(nBas,nspin),Z(nBas,nspin), & - eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), & + allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & + Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & + Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & + SigT(nBas,nspin),Z(nBas,nspin),eGT(nBas,nspin), & + eOld(nBas,nspin),error_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 iblock = 3 -! iblock = 1 ! 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 !---------------------------------------------- - EcGM = 0d0 - SigT(:,:) = 0d0 - Z(:,:) = 0d0 - !alpha-beta block iblock = 3 @@ -199,13 +194,7 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & 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,& - rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) - - 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(:,:)) + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z) !---------------------------------------------- ! Solve the quasi-particle equation diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index f69054b..4595bf5 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -248,10 +248,6 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute T-matrix version of the self-energy !---------------------------------------------- - EcGM = 0d0 - SigT(:,:,:) = 0d0 - Z(:,:) = 0d0 - !alpha-beta block iblock = 3 @@ -273,13 +269,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & rho1bb,X2bb,Y2bb,rho2bb) 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) - - 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(:,:)) + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT,Z) ! 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 F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) & - + SigTp(:,:,ispin) + + SigTp(:,:,ispin) end do ! Compute commutator and convergence criteria