From d6eb90df47a9f1aebf0d89566d0e7a693027c454 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 12 Jul 2023 14:13:45 +0200 Subject: [PATCH] moving calculations of Z in Sig --- src/GT/G0T0eh.f90 | 15 ++-- src/GT/G0T0pp.f90 | 20 +++-- src/GT/GTeh_renormalization_factor.f90 | 61 -------------- src/GT/GTeh_self_energy.f90 | 39 ++++++--- src/GT/GTeh_self_energy_diag.f90 | 33 ++++++-- src/GT/GTpp_renormalization_factor.f90 | 51 ------------ src/GT/GTpp_self_energy.f90 | 43 ++++++---- src/GT/GTpp_self_energy_diag.f90 | 41 +++++++--- src/GT/evGTeh.f90 | 15 ++-- src/GT/evGTpp.f90 | 28 +++---- src/GT/qsGTeh.f90 | 27 +++---- src/GT/qsGTpp.f90 | 36 ++++----- src/GW/G0W0.f90 | 3 +- src/GW/GW_renormalization_factor.f90 | 72 ----------------- src/GW/GW_self_energy.f90 | 105 +++++++++++++++---------- src/GW/GW_self_energy_diag.f90 | 39 ++++++--- src/GW/evGW.f90 | 3 +- src/GW/qsGW.f90 | 3 +- 18 files changed, 263 insertions(+), 371 deletions(-) delete mode 100644 src/GT/GTeh_renormalization_factor.f90 delete mode 100644 src/GT/GTpp_renormalization_factor.f90 delete mode 100644 src/GW/GW_renormalization_factor.f90 diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 75d14db..fb37dfe 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -53,7 +53,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD double precision :: EcppBSE(nspin) double precision :: EcGM double precision,allocatable :: SigX(:) - double precision,allocatable :: SigC(:) + double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) @@ -99,7 +99,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Memory allocation - allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + allocate(Sig(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas)) !-------------------! @@ -124,13 +124,12 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD if(regularize) then -! call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) +! call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,Sig) ! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else - call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) - call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,Z) + call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z) end if @@ -138,7 +137,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Solve the quasi-particle equation ! !-----------------------------------! - eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) + eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) ! Linearized or graphical solution? @@ -170,11 +169,11 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Dump results ! !--------------! - call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM) ! Deallocate memory -! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) +! deallocate(Sig,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) ! Plot stuff diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index db9589f..c0f29fa 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -60,7 +60,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp double precision,allocatable :: Y2ab(:,:),Y2aa(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:) double precision,allocatable :: SigX(:) - double precision,allocatable :: SigT(:) + double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: eGTlin(:) @@ -92,7 +92,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), & Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & - SigX(nBas),SigT(nBas),Z(nBas),eGTlin(nBas)) + SigX(nBas),Sig(nBas),Z(nBas),eGTlin(nBas)) !---------------------------------------------- ! alpha-beta block @@ -127,7 +127,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp !---------------------------------------------- EcGM = 0d0 - SigT(:) = 0d0 + Sig(:) = 0d0 Z(:) = 0d0 iblock = 3 @@ -136,13 +136,12 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp if(regularize) then - call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,SigT) + call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,Sig) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,Z) else - call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,SigT) - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,Z) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOab,nVVab,eHF,Om1ab,rho1ab,Om2ab,rho2ab,EcGM,Sig,Z) end if @@ -152,13 +151,12 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp if(regularize) then - call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,SigT) + call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,Sig) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,Z) else - call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,SigT) - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,Z) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,EcGM,Sig,Z) end if @@ -174,7 +172,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp ! Solve the quasi-particle equation !---------------------------------------------- - eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + SigT(:) - Vxc(:)) + eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) if(linearize) then @@ -210,7 +208,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) - call print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) + call print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) ! Perform BSE calculation diff --git a/src/GT/GTeh_renormalization_factor.f90 b/src/GT/GTeh_renormalization_factor.f90 deleted file mode 100644 index 75ad573..0000000 --- a/src/GT/GTeh_renormalization_factor.f90 +++ /dev/null @@ -1,61 +0,0 @@ -subroutine GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,Z) - -! Compute renormalization factor for GW - - 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 - integer,intent(in) :: nS - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rhoL(nBas,nBas,nS) - double precision,intent(in) :: rhoR(nBas,nBas,nS) - -! Local variables - - integer :: p,i,a,m - double precision :: eps - -! Output variables - - double precision,intent(out) :: Z(nBas) - -! Initialize - - Z(:) = 0d0 - -! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do m=1,nS - eps = e(p) - e(i) + Om(m) - Z(p) = Z(p) - rhoL(i,p,m)*rhoR(i,p,m)*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - -! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do m=1,nS - eps = e(p) - e(a) - Om(m) - Z(p) = Z(p) - rhoL(p,a,m)*rhoR(p,a,m)*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - -! Compute renormalization factor from derivative of SigC - - Z(:) = 1d0/(1d0 - Z(:)) - -end subroutine diff --git a/src/GT/GTeh_self_energy.f90 b/src/GT/GTeh_self_energy.f90 index af06e80..4a7fd43 100644 --- a/src/GT/GTeh_self_energy.f90 +++ b/src/GT/GTeh_self_energy.f90 @@ -1,6 +1,6 @@ -subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) +subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig,Z) -! Compute correlation part of the self-energy for GTeh +! Compute correlation part of the self-energy for GTeh and the renormalization factor implicit none include 'parameters.h' @@ -24,16 +24,18 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) integer :: i,j,a,b integer :: p,q,r integer :: m - double precision :: eps + double precision :: num,eps ! Output variables double precision,intent(out) :: EcGM - double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: Sig(nBas,nBas) + double precision,intent(out) :: Z(nBas) ! Initialize - SigC(:,:) = 0d0 + Sig(:,:) = 0d0 + Z(:) = 0d0 !----------------! ! GW self-energy ! @@ -42,16 +44,20 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) ! Occupied part of the correlation self-energy !$OMP PARALLEL & -!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & -!$OMP PRIVATE(m,i,q,p,eps) & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) & +!$OMP PRIVATE(m,i,q,p,num,eps) & !$OMP DEFAULT(NONE) !$OMP DO do q=nC+1,nBas-nR do p=nC+1,nBas-nR do m=1,nS do i=nC+1,nO + eps = e(p) - e(i) + Om(m) - SigC(p,q) = SigC(p,q) + rhoL(i,p,m)*rhoR(i,q,m)*eps/(eps**2 + eta**2) + num = rhoL(i,p,m)*rhoR(i,q,m) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -62,16 +68,20 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) ! Virtual part of the correlation self-energy !$OMP PARALLEL & -!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & -!$OMP PRIVATE(m,a,q,p,eps) & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) & +!$OMP PRIVATE(m,a,q,p,num,eps) & !$OMP DEFAULT(NONE) !$OMP DO do q=nC+1,nBas-nR do p=nC+1,nBas-nR do m=1,nS do a=nO+1,nBas-nR + eps = e(p) - e(a) - Om(m) - SigC(p,q) = SigC(p,q) + rhoL(p,a,m)*rhoR(q,a,m)*eps/(eps**2 + eta**2) + num = rhoL(p,a,m)*rhoR(q,a,m) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -86,9 +96,14 @@ subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) do a=nO+1,nBas-nR do i=nC+1,nO eps = e(a) - e(i) + Om(m) - EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2) + num = rhoL(i,a,m)*rhoR(i,a,m) + EcGM = EcGM - num*eps/(eps**2 + eta**2) end do end do end do +! Compute renormalization factor from derivative + + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 index 4dd666f..f48c06c 100644 --- a/src/GT/GTeh_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -1,6 +1,6 @@ -subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) +subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig,Z) -! Compute diagonal of the correlation part of the self-energy +! Compute diagonal of the correlation part of the self-energy and the renormalization factor implicit none include 'parameters.h' @@ -22,16 +22,18 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig ! Local variables integer :: i,a,p,q,m - double precision :: eps + double precision :: num,eps ! Output variables - double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: Sig(nBas) + double precision,intent(out) :: Z(nBas) double precision,intent(out) :: EcGM ! Initialize - SigC(:) = 0d0 + Sig(:) = 0d0 + Z(:) = 0d0 !----------------------------- ! GW self-energy @@ -42,8 +44,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do p=nC+1,nBas-nR do i=nC+1,nO do m=1,nS + eps = e(p) - e(i) + Om(m) - SigC(p) = SigC(p) + rhoL(i,p,m)*rhoR(i,p,m)*eps/(eps**2 + eta**2) + num = rhoL(i,p,m)*rhoR(i,p,m) + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -53,8 +59,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do p=nC+1,nBas-nR do a=nO+1,nBas-nR do m=1,nS + eps = e(p) - e(a) - Om(m) - SigC(p) = SigC(p) + rhoL(p,a,m)*rhoR(p,a,m)*eps/(eps**2 + eta**2) + num = rhoL(p,a,m)*rhoR(p,a,m) + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -65,10 +75,17 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do i=nC+1,nO do a=nO+1,nBas-nR do m=1,nS + eps = e(a) - e(i) + Om(m) - EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2) + num = rhoL(i,a,m)*rhoR(i,a,m) + EcGM = EcGM - num*eps/(eps**2 + eta**2) + end do end do end do +! Compute renormalization factor from derivative + + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GT/GTpp_renormalization_factor.f90 b/src/GT/GTpp_renormalization_factor.f90 deleted file mode 100644 index cea316a..0000000 --- a/src/GT/GTpp_renormalization_factor.f90 +++ /dev/null @@ -1,51 +0,0 @@ -subroutine GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,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,nO,nV,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,p,cd,kl - double precision :: eps - -! Output variables - - double precision,intent(out) :: Z(nBas) - -! Occupied part of the T-matrix self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do cd=1,nVV - eps = e(p) + e(i) - Omega1(cd) - Z(p) = Z(p) - rho1(p,i,cd)**2*(eps/(eps**2 + eta**2))**2 - enddo - enddo - enddo - -! Virtual part of the T-matrix self-energy - - do p=nC+1,nBas-nR - do a=1,nV-nR - do kl=1,nOO - eps = e(p) + e(nO+a) - Omega2(kl) - Z(p) = Z(p) - rho2(p,nO+a,kl)**2*(eps/(eps**2 + eta**2))**2 - enddo - enddo - enddo - -end subroutine diff --git a/src/GT/GTpp_self_energy.f90 b/src/GT/GTpp_self_energy.f90 index f8cd2ad..27113e7 100644 --- a/src/GT/GTpp_self_energy.f90 +++ b/src/GT/GTpp_self_energy.f90 @@ -1,6 +1,6 @@ -subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) +subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) -! Compute the correlation part of the T-matrix self-energy +! Compute the correlation part of the T-matrix self-energy and the renormalization factor implicit none include 'parameters.h' @@ -16,20 +16,21 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: rho2(nBas,nBas,nOO) ! 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 - double precision,intent(inout):: SigT(nBas,nBas) + double precision,intent(inout):: Sig(nBas,nBas) + double precision,intent(inout):: Z(nBas) !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -39,8 +40,12 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh do q=nC+1,nBas-nR do i=nC+1,nO do cd=1,nVV - eps = e(p) + e(i) - Omega1(cd) - SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2) + + eps = e(p) + e(i) - Om1(cd) + num = rho1(p,i,cd)*rho1(q,i,cd) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo enddo enddo @@ -54,8 +59,12 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh do q=nC+1,nBas-nR do a=nO+1,nBas-nR do kl=1,nOO - eps = e(p) + e(a) - Omega2(kl) - SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2) + + eps = e(p) + e(a) - Om2(kl) + num = rho2(p,a,kl)*rho2(q,a,kl) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo enddo enddo @@ -68,8 +77,11 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh do i=nC+1,nO do j=nC+1,nO do cd=1,nVV - eps = e(i) + e(j) - Omega1(cd) - EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + + eps = e(i) + e(j) - Om1(cd) + num = rho1(i,j,cd)*rho1(i,j,cd) + EcGM = EcGM + num*eps/(eps**2 + eta**2) + enddo enddo enddo @@ -77,8 +89,11 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rh do a=nO+1,nBas-nR do b=nO+1,nBas-nR do kl=1,nOO - eps = e(a) + e(b) - Omega2(kl) - EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + + eps = e(a) + e(b) - Om2(kl) + num = rho2(a,b,kl)*rho2(a,b,kl) + EcGM = EcGM - num*eps/(eps**2 + eta**2) + enddo enddo enddo diff --git a/src/GT/GTpp_self_energy_diag.f90 b/src/GT/GTpp_self_energy_diag.f90 index e797b44..6d9e2ed 100644 --- a/src/GT/GTpp_self_energy_diag.f90 +++ b/src/GT/GTpp_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) +subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -16,20 +16,21 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: rho2(nBas,nBas,nOO) ! Local variables integer :: i,j,a,b,p,cd,kl - double precision :: eps + double precision :: num,eps ! Output variables double precision,intent(inout) :: EcGM - double precision,intent(inout) :: SigT(nBas) + double precision,intent(inout) :: Sig(nBas) + double precision,intent(inout) :: Z(nBas) !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -38,8 +39,12 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg do p=nC+1,nBas-nR do i=nC+1,nO do cd=1,nVV - eps = e(p) + e(i) - Omega1(cd) - SigT(p) = SigT(p) + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + + eps = e(p) + e(i) - Om1(cd) + num = rho1(p,i,cd)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo enddo enddo @@ -51,8 +56,12 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg do p=nC+1,nBas-nR do a=nO+1,nBas-nR do kl=1,nOO - eps = e(p) + e(a) - Omega2(kl) - SigT(p) = SigT(p) + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + + eps = e(p) + e(a) - Om2(kl) + num = rho2(p,a,kl)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo enddo enddo @@ -64,8 +73,11 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg do i=nC+1,nO do j=nC+1,nO do cd=1,nVV - eps = e(i) + e(j) - Omega1(cd) - EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + + eps = e(i) + e(j) - Om1(cd) + num = rho1(i,j,cd)**2 + EcGM = EcGM + num*eps/(eps**2 + eta**2) + enddo enddo enddo @@ -73,8 +85,11 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omeg do a=nO+1,nBas-nR do b=nO+1,nBas-nR do kl=1,nOO - eps = e(a) + e(b) - Omega2(kl) - EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + + eps = e(a) + e(b) - Om2(kl) + num = rho2(a,b,kl)**2 + EcGM = EcGM - num*eps/(eps**2 + eta**2) + enddo enddo enddo diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index 6fad6f3..07ee8b1 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -68,7 +68,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) double precision,allocatable :: SigX(:) - double precision,allocatable :: SigC(:) + double precision,allocatable :: Sig(:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) @@ -106,7 +106,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Memory allocation - allocate(eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + allocate(eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),Sig(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas)) ! Compute the exchange part of the self-energy @@ -146,19 +146,18 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, if(regularize) then -! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) +! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,Sig) ! call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) else - call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) - call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) + call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z) end if ! Solve the quasi-particle equation - eGTlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + eGTlin(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) ! Linearized or graphical solution? @@ -184,7 +183,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Print results - call print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + call print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM) ! Linear mixing or DIIS extrapolation @@ -232,7 +231,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Deallocate memory - deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error_diis,e_diis) + deallocate(eOld,Z,Sig,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error_diis,e_diis) ! Perform BSE calculation diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index 2167cf5..e6a6cae 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -71,7 +71,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) double precision,allocatable :: SigX(:) - double precision,allocatable :: SigT(:) + double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) ! Output variables @@ -100,7 +100,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & - eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigT(nBas), & + eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),Sig(nBas), & error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Compute the exchange part of the self-energy @@ -155,30 +155,20 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & !---------------------------------------------- EcGM = 0d0 - SigT(:) = 0d0 + Sig(:) = 0d0 Z(:) = 0d0 iblock = 3 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, & - X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Om1s,rho1s,Om2s,rho2s,EcGM,SigT) - - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Om1s,rho1s,Om2s,rho2s,Z) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig,Z) iblock = 4 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, & - X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Om1t,rho1t,Om2t,rho2t,EcGM,SigT) - - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Om1t,rho1t,Om2t,rho2t,Z) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) Z(:) = 1d0/(1d0 - Z(:)) @@ -188,7 +178,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Solve the quasi-particle equation !---------------------------------------------- - eGT(:) = eHF(:) + SigX(:) + SigT(:) - Vxc(:) + eGT(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) ! Convergence criteria @@ -198,7 +188,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Dump results !---------------------------------------------- - call print_evGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) + call print_evGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) ! DIIS extrapolation diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsGTeh.f90 index e0b6161..343329d 100644 --- a/src/GT/qsGTeh.f90 +++ b/src/GT/qsGTeh.f90 @@ -89,9 +89,9 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, double precision,allocatable :: Fp(:,:) double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) - double precision,allocatable :: SigC(:,:) - double precision,allocatable :: SigCp(:,:) - double precision,allocatable :: SigCm(:,:) + double precision,allocatable :: Sig(:,:) + double precision,allocatable :: Sigp(:,:) + double precision,allocatable :: Sigm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -129,7 +129,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Memory allocation allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & + J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), & OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) @@ -180,26 +180,25 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, if(regularize) then -! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) +! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig) ! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) else - call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) - call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) + call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,Sig,Z) endif ! Make correlation self-energy Hermitian and transform it back to AO basis - SigCp = 0.5d0*(SigC + transpose(SigC)) - SigCm = 0.5d0*(SigC - transpose(SigC)) + Sigp = 0.5d0*(Sig + transpose(Sig)) + Sigm = 0.5d0*(Sig - transpose(Sig)) - call MOtoAO_transform(nBas,S,c,SigCp) + call MOtoAO_transform(nBas,S,c,Sigp) ! Solve the quasi-particle equation - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:) ! Compute commutator and convergence criteria @@ -220,7 +219,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGT) c = matmul(X,cp) - SigCp = matmul(transpose(c),matmul(SigCp,c)) + Sigp = matmul(transpose(c),matmul(Sigp,c)) ! Compute new density matrix in the AO basis @@ -258,7 +257,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) enddo !------------------------------------------------------------------------ @@ -281,7 +280,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 8a1a1a8..776e8a8 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -89,9 +89,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T double precision,allocatable :: Fp(:,:) double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) - double precision,allocatable :: SigT(:,:) - double precision,allocatable :: SigTp(:,:) - double precision,allocatable :: SigTm(:,:) + double precision,allocatable :: Sig(:,:) + double precision,allocatable :: Sigp(:,:) + double precision,allocatable :: Sigm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -137,7 +137,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Memory allocation allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigT(nBas,nBas),SigTp(nBas,nBas),SigTm(nBas,nBas),Z(nBas), & + J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & @@ -201,7 +201,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Compute correlation part of the self-energy EcGM = 0d0 - SigT(:,:) = 0d0 + Sig(:,:) = 0d0 Z(:) = 0d0 iblock = 3 @@ -210,15 +210,13 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T if(regularize) then - call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,SigT) + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,Z) else - call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,SigT) - - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,Z) + call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s,EcGM,Sig,Z) end if @@ -228,16 +226,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T if(regularize) then - call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,SigT) + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig) call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & Om1t,rho1t,Om2t,rho2t,Z) else - call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,SigT) - - call GTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,Z) + call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) end if @@ -245,14 +241,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Make correlation self-energy Hermitian and transform it back to AO basis - SigTp = 0.5d0*(SigT + transpose(SigT)) - SigTm = 0.5d0*(SigT - transpose(SigT)) + Sigp = 0.5d0*(Sig + transpose(Sig)) + Sigm = 0.5d0*(Sig - transpose(Sig)) - call MOtoAO_transform(nBas,S,c,SigTp) + call MOtoAO_transform(nBas,S,c,Sigp) ! Solve the quasi-particle equation - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigTp(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:) ! Compute commutator and convergence criteria @@ -273,7 +269,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGT) c = matmul(X,cp) - SigTp = matmul(transpose(c),matmul(SigTp,c)) + Sigp = matmul(transpose(c),matmul(Sigp,c)) ! Compute new density matrix in the AO basis @@ -311,7 +307,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigTp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call print_qsGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) enddo !------------------------------------------------------------------------ @@ -334,7 +330,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 99addb8..70ec2a9 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -144,8 +144,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD else - call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC,Z) end if diff --git a/src/GW/GW_renormalization_factor.f90 b/src/GW/GW_renormalization_factor.f90 deleted file mode 100644 index 3820235..0000000 --- a/src/GW/GW_renormalization_factor.f90 +++ /dev/null @@ -1,72 +0,0 @@ -subroutine GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) - -! Compute renormalization factor for GW - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: COHSEX - 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) :: nS - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: p,i,a,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: Z(nBas) - -! Initialize - - Z(:) = 0d0 - -! static COHSEX approximation - - if(COHSEX) then - - Z(:) = 1d0 - return - - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - end do - end do - - end if - -! Compute renormalization factor from derivative of SigC - - Z(:) = 1d0/(1d0 - Z(:)) - -end subroutine diff --git a/src/GW/GW_self_energy.f90 b/src/GW/GW_self_energy.f90 index 91917e2..687385f 100644 --- a/src/GW/GW_self_energy.f90 +++ b/src/GW/GW_self_energy.f90 @@ -1,6 +1,6 @@ -subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) +subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) -! Compute correlation part of the self-energy +! Compute correlation part of the self-energy and the renormalization factor implicit none include 'parameters.h' @@ -24,16 +24,18 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) integer :: i,j,a,b integer :: p,q,r integer :: jb - double precision :: eps + double precision :: num,eps ! Output variables double precision,intent(out) :: EcGM - double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: Sig(nBas,nBas) + double precision,intent(out) :: Z(nBas) ! Initialize - SigC(:,:) = 0d0 + Sig(:,:) = 0d0 + Z(:) = 0d0 !-----------------------------! ! COHSEX static approximation ! @@ -47,7 +49,7 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) do q=nC+1,nBas-nR do i=nC+1,nO do jb=1,nS - SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) + Sig(p,q) = Sig(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) end do end do end do @@ -59,7 +61,7 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) do q=nC+1,nBas-nR do r=nC+1,nBas-nR do jb=1,nS - SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) + Sig(p,q) = Sig(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) end do end do end do @@ -67,9 +69,11 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM + 0.5d0*SigC(i,i) + EcGM = EcGM + 0.5d0*Sig(i,i) end do + Z(:) = 0d0 + else !----------------! @@ -78,43 +82,51 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Occupied part of the correlation self-energy -!$OMP PARALLEL & -!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & -!$OMP PRIVATE(jb,i,q,p,eps) & -!$OMP DEFAULT(NONE) -!$OMP DO -do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do jb=1,nS - do i=nC+1,nO - eps = e(p) - e(i) + Omega(jb) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2) - end do - end do - end do -end do -!$OMP END DO -!$OMP END PARALLEL + !$OMP PARALLEL & + !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & + !$OMP PRIVATE(jb,i,q,p,eps,num) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do jb=1,nS + do i=nC+1,nO + + eps = e(p) - e(i) + Omega(jb) + num = 2d0*rho(p,i,jb)*rho(q,i,jb) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL ! Virtual part of the correlation self-energy -!$OMP PARALLEL & -!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & -!$OMP PRIVATE(jb,a,q,p,eps) & -!$OMP DEFAULT(NONE) -!$OMP DO -do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do jb=1,nS - do a=nO+1,nBas-nR - eps = e(p) - e(a) - Omega(jb) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2) - end do - end do - end do -end do -!$OMP END DO -!$OMP END PARALLEL + !$OMP PARALLEL & + !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & + !$OMP PRIVATE(jb,a,q,p,eps,num) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do jb=1,nS + do a=nO+1,nBas-nR + + eps = e(p) - e(a) - Omega(jb) + num = 2d0*rho(p,a,jb)*rho(q,a,jb) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL ! Galitskii-Migdal correlation energy @@ -122,12 +134,19 @@ end do do jb=1,nS do a=nO+1,nBas-nR do i=nC+1,nO + eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) + num = 4d0*rho(a,i,jb)*rho(a,i,jb) + EcGM = EcGM - num*eps/(eps**2 + eta**2) + end do end do end do end if +! Compute renormalization factor from derivative + + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GW/GW_self_energy_diag.f90 b/src/GW/GW_self_energy_diag.f90 index e02f5fa..735591e 100644 --- a/src/GW/GW_self_energy_diag.f90 +++ b/src/GW/GW_self_energy_diag.f90 @@ -1,6 +1,6 @@ -subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) +subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) -! Compute diagonal of the correlation part of the self-energy +! Compute diagonal of the correlation part of the self-energy and the renormalization factor implicit none include 'parameters.h' @@ -22,16 +22,18 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S ! Local variables integer :: i,a,p,q,jb - double precision :: eps + double precision :: num,eps ! Output variables - double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: Sig(nBas) + double precision,intent(out) :: Z(nBas) double precision,intent(out) :: EcGM ! Initialize - SigC(:) = 0d0 + Sig(:) = 0d0 + Z(:) = 0d0 !----------------------------- ! COHSEX static self-energy @@ -44,7 +46,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S do p=nC+1,nBas-nR do i=nC+1,nO do jb=1,nS - SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/Omega(jb) + Sig(p) = Sig(p) + 4d0*rho(p,i,jb)**2/Omega(jb) end do end do end do @@ -54,7 +56,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S do p=nC+1,nBas-nR do q=nC+1,nBas-nR do jb=1,nS - SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb) + Sig(p) = Sig(p) - 2d0*rho(p,q,jb)**2/Omega(jb) end do end do end do @@ -63,7 +65,7 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM - SigC(i) + EcGM = EcGM - Sig(i) end do !----------------------------- @@ -77,8 +79,12 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S do p=nC+1,nBas-nR do i=nC+1,nO do jb=1,nS + eps = e(p) - e(i) + Omega(jb) - SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) + num = 2d0*rho(p,i,jb)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -88,8 +94,12 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S do p=nC+1,nBas-nR do a=nO+1,nBas-nR do jb=1,nS + eps = e(p) - e(a) - Omega(jb) - SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) + num = 2d0*rho(p,a,jb)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do end do end do @@ -100,12 +110,19 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S do i=nC+1,nO do a=nO+1,nBas-nR do jb=1,nS + eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) + num = 4d0*rho(a,i,jb)**2 + EcGM = EcGM - num*eps/(eps**2 + eta**2) + end do end do end do end if +! Compute renormalization factor from derivative + + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 7a518ca..28f7b5f 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -167,8 +167,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, else - call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z) end if diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index cc645bc..106709b 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -192,8 +192,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, else - call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call GW_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z) endif