From 4fec6e88f5b8132e9dd0b67a9dbe2289135df1d7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 27 Jul 2023 19:17:20 +0200 Subject: [PATCH] clean up GTpp --- input/methods | 2 +- input/options | 2 +- src/GT/G0T0pp.f90 | 142 +++++++++++++------------------ src/GT/GT.f90 | 4 +- src/GT/GTpp_QP_graph.f90 | 19 +++-- src/GT/GTpp_SigC.f90 | 41 ++++++--- src/GT/GTpp_dSigC.f90 | 42 ++++++--- src/GT/GTpp_self_energy.f90 | 80 ++++++++++++----- src/GT/GTpp_self_energy_diag.f90 | 84 ++++++++++++------ src/GT/evGTpp.f90 | 45 +++++----- src/GT/qsGTpp.f90 | 46 ++-------- 11 files changed, 278 insertions(+), 229 deletions(-) diff --git a/input/methods b/input/methods index 4ef876d..bb2e24b 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F T F F + T T F F F F # * unrestricted version available diff --git a/input/options b/input/options index 0b1a49b..dd8f88a 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg 256 0.00001 T 5 T 0.0 F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 256 0.00001 T 5 T 0.0 F F + 256 0.00001 T 5 F 0.0 F F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index d4bb341..92082ee 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -41,22 +41,22 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp integer :: ispin integer :: iblock - integer :: nOOab,nOOaa - integer :: nVVab,nVVaa + integer :: nOOs,nOOt + integer :: nVVs,nVVt double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcGM double precision,allocatable :: Bpp(:,:) double precision,allocatable :: Cpp(:,:) double precision,allocatable :: Dpp(:,:) - double precision,allocatable :: Om1ab(:),Om1aa(:) - double precision,allocatable :: X1ab(:,:),X1aa(:,:) - double precision,allocatable :: Y1ab(:,:),Y1aa(:,:) - double precision,allocatable :: rho1ab(:,:,:),rho1aa(:,:,:) - double precision,allocatable :: Om2ab(:),Om2aa(:) - double precision,allocatable :: X2ab(:,:),X2aa(:,:) - double precision,allocatable :: Y2ab(:,:),Y2aa(:,:) - double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:) + double precision,allocatable :: Om1s(:),Om1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Om2s(:),Om2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: eGT(:) @@ -74,20 +74,20 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ! Dimensions of the pp-RPA linear reponse matrices - nOOab = nO*nO - nVVab = nV*nV + nOOs = nO*nO + nVVs = nV*nV - nOOaa = nO*(nO - 1)/2 - nVVaa = nV*(nV - 1)/2 + nOOt = nO*(nO - 1)/2 + nVVt = nV*(nV - 1)/2 ! Memory allocation - allocate(Om1ab(nVVab),X1ab(nVVab,nVVab),Y1ab(nOOab,nVVab), & - Om2ab(nOOab),X2ab(nVVab,nOOab),Y2ab(nOOab,nOOab), & - rho1ab(nBas,nBas,nVVab),rho2ab(nBas,nBas,nOOab), & - Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), & - Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & - rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas)) !---------------------------------------------- @@ -99,18 +99,18 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ! Compute linear response - allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eHF,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp) - call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) - call print_excitation('pp-RPA (N+2)',iblock,nVVab,Om1ab(:)) - call print_excitation('pp-RPA (N-2)',iblock,nOOab,Om2ab(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVs,Om1s(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOs,Om2s(:)) !---------------------------------------------- ! alpha-alpha block @@ -121,78 +121,50 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ! Compute linear response - allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eHF,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp) - call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) - call print_excitation('pp-RPA (N+2)',iblock,nVVaa,Om1aa(:)) - call print_excitation('pp-RPA (N-2)',iblock,nOOaa,Om2aa(:)) + call print_excitation('ppRPA (N+2) ',iblock,nVVt,Om1t) + call print_excitation('ppRPA (N-2) ',iblock,nOOt,Om2t) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- - EcGM = 0d0 - Sig(:) = 0d0 - Z(:) = 0d0 - iblock = 3 - - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,ERI_MO,X1ab,Y1ab,rho1ab,X2ab,Y2ab,rho2ab) - - if(regularize) then - - 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,Sig,Z) - - end if + call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) 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,nOOaa,nVVaa,ERI_MO,X1aa,Y1aa,rho1aa,X2aa,Y2aa,rho2aa) - - if(regularize) then - - 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,Sig,Z) - - end if - - Z(:) = 1d0/(1d0 - Z(:)) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) !---------------------------------------------- ! Solve the quasi-particle equation !---------------------------------------------- - - eGTlin(:) = eHF(:) + Z(:)*Sig(:) if(linearize) then write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' write(*,*) - eGT(:) = eGTlin(:) + eGT(:) = eHF(:) + Z(:)*Sig(:) else write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT) + call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,eHF,eGT) end if @@ -205,26 +177,26 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp ispin = 1 iblock = 3 - allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) - call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) ispin = 2 iblock = 4 - allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) - call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) @@ -237,8 +209,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp if(dophBSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & - Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int,eHF,eGT,EcBSE) if(exchange_kernel) then @@ -274,8 +246,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp end if call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - nOOab,nVVab,nOOaa,nVVaa,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa, & - Om2aa,X2aa,Y2aa,rho1aa,rho2aa,ERI_MO,eHF,eGT,EcBSE) + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eHF,eGT,EcBSE) if(exchange_kernel) then @@ -299,8 +271,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp if(doppBSE) then - call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & - Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, & + Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int,eHF,eGT,EcBSE) write(*,*) diff --git a/src/GT/GT.f90 b/src/GT/GT.f90 index db7fdca..702ebc8 100644 --- a/src/GT/GT.f90 +++ b/src/GT/GT.f90 @@ -108,11 +108,11 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted call wall_time(start_GT) if(unrestricted) then call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & PHF,cHF,epsHF) else call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,PHF,cHF,epsHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,PHF,cHF,epsHF) end if call wall_time(end_GT) diff --git a/src/GT/GTpp_QP_graph.f90 b/src/GT/GTpp_QP_graph.f90 index d79291a..1e0949b 100644 --- a/src/GT/GTpp_QP_graph.f90 +++ b/src/GT/GTpp_QP_graph.f90 @@ -1,4 +1,5 @@ -subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eGT) +subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,eGTlin,eGT) implicit none include 'parameters.h' @@ -9,15 +10,15 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTl integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: eta double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) double precision,intent(in) :: eGTlin(nBas) @@ -54,8 +55,8 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTl nIt = nIt + 1 - sigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2) - dsigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2) + sigC = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) + dsigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t) write (*,*) sigC f = w - eHF(p) - sigC df = 1d0 - dsigC diff --git a/src/GT/GTpp_SigC.f90 b/src/GT/GTpp_SigC.f90 index c151d65..c511524 100644 --- a/src/GT/GTpp_SigC.f90 +++ b/src/GT/GTpp_SigC.f90 @@ -1,4 +1,5 @@ -double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2) +double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t) ! Compute diagonal of the correlation part of the self-energy @@ -15,12 +16,14 @@ double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1, integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO,nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt + double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -36,10 +39,17 @@ double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1, !---------------------------------------------- do i=nC+1,nO - do cd=1,nVV - eps = w + e(i) - Om1(cd) - GTpp_SigC = GTpp_SigC + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + + do cd=1,nVVs + eps = w + e(i) - Om1s(cd) + GTpp_SigC = GTpp_SigC + rho1s(p,i,cd)**2*eps/(eps**2 + eta**2) enddo + + do cd=1,nVVt + eps = w + e(i) - Om1t(cd) + GTpp_SigC = GTpp_SigC + rho1t(p,i,cd)**2*eps/(eps**2 + eta**2) + enddo + enddo !---------------------------------------------- @@ -47,10 +57,17 @@ double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1, !---------------------------------------------- do a=nO+1,nBas-nR - do kl=1,nOO - eps = w + e(a) - Om2(kl) - GTpp_SigC = GTpp_SigC + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + + do kl=1,nOOs + eps = w + e(a) - Om2s(kl) + GTpp_SigC = GTpp_SigC + rho2s(p,a,kl)**2*eps/(eps**2 + eta**2) enddo + + do kl=1,nOOt + eps = w + e(a) - Om2t(kl) + GTpp_SigC = GTpp_SigC + rho2t(p,a,kl)**2*eps/(eps**2 + eta**2) + enddo + enddo end function diff --git a/src/GT/GTpp_dSigC.f90 b/src/GT/GTpp_dSigC.f90 index cced4dd..1251ab4 100644 --- a/src/GT/GTpp_dSigC.f90 +++ b/src/GT/GTpp_dSigC.f90 @@ -1,4 +1,5 @@ -double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2) +double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t) ! Compute diagonal of the correlation part of the self-energy @@ -15,13 +16,14 @@ double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1 integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt + double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -37,12 +39,17 @@ double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1 !---------------------------------------------- do i=nC+1,nO - do cd=1,nVV - - eps = w + e(i) - Om1(cd) - GTpp_dSigC= GTpp_dSigC - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + do cd=1,nVVs + eps = w + e(i) - Om1s(cd) + GTpp_dSigC = GTpp_dSigC - rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo + + do cd=1,nVVt + eps = w + e(i) - Om1t(cd) + GTpp_dSigC = GTpp_dSigC - rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo !---------------------------------------------- @@ -51,10 +58,17 @@ double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1 !---------------------------------------------- do a=nO+1,nBas-nR - do kl=1,nOO - eps = w + e(a) - Om2(kl) - GTpp_dSigC = GTpp_dSigC - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + do kl=1,nOOs + eps = w + e(a) - Om2s(kl) + GTpp_dSigC = GTpp_dSigC - rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo + + do kl=1,nOOt + eps = w + e(a) - Om2t(kl) + GTpp_dSigC = GTpp_dSigC - rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo end function diff --git a/src/GT/GTpp_self_energy.f90 b/src/GT/GTpp_self_energy.f90 index 27113e7..49cdefb 100644 --- a/src/GT/GTpp_self_energy.f90 +++ b/src/GT/GTpp_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) +subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Compute the correlation part of the T-matrix self-energy and the renormalization factor @@ -13,13 +13,13 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcG integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -32,6 +32,12 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcG double precision,intent(inout):: Sig(nBas,nBas) double precision,intent(inout):: Z(nBas) +! Initialization + + Sig(:,:) = 0d0 + Z(:) = 0d0 + EcGM = 0d0 + !---------------------------------------------- ! Occupied part of the T-matrix self-energy !---------------------------------------------- @@ -39,14 +45,21 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcG do p=nC+1,nBas-nR do q=nC+1,nBas-nR do i=nC+1,nO - do cd=1,nVV - eps = e(p) + e(i) - Om1(cd) - num = rho1(p,i,cd)*rho1(q,i,cd) + do cd=1,nVVs + eps = e(p) + e(i) - Om1s(cd) + num = rho1s(p,i,cd)*rho1s(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 + + do cd=1,nVVt + eps = e(p) + e(i) - Om1t(cd) + num = rho1t(p,i,cd)*rho1t(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 enddo @@ -58,14 +71,21 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcG do p=nC+1,nBas-nR do q=nC+1,nBas-nR do a=nO+1,nBas-nR - do kl=1,nOO - eps = e(p) + e(a) - Om2(kl) - num = rho2(p,a,kl)*rho2(q,a,kl) + do kl=1,nOOs + eps = e(p) + e(a) - Om2s(kl) + num = rho2s(p,a,kl)*rho2s(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 + + do kl=1,nOOt + eps = e(p) + e(a) - Om2t(kl) + num = rho2t(p,a,kl)*rho2t(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 enddo @@ -76,26 +96,40 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcG do i=nC+1,nO do j=nC+1,nO - do cd=1,nVV - eps = e(i) + e(j) - Om1(cd) - num = rho1(i,j,cd)*rho1(i,j,cd) + do cd=1,nVVs + eps = e(i) + e(j) - Om1s(cd) + num = rho1s(i,j,cd)*rho1s(i,j,cd) EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + + do cd=1,nVVt + eps = e(i) + e(j) - Om1t(cd) + num = rho1t(i,j,cd)*rho1t(i,j,cd) + EcGM = EcGM + num*eps/(eps**2 + eta**2) + enddo + enddo enddo do a=nO+1,nBas-nR do b=nO+1,nBas-nR - do kl=1,nOO - eps = e(a) + e(b) - Om2(kl) - num = rho2(a,b,kl)*rho2(a,b,kl) + do kl=1,nOOs + eps = e(a) + e(b) - Om2s(kl) + num = rho2s(a,b,kl)*rho2s(a,b,kl) EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + + do kl=1,nOOt + eps = e(a) + e(b) - Om2t(kl) + num = rho2t(a,b,kl)*rho2t(a,b,kl) + EcGM = EcGM - num*eps/(eps**2 + eta**2) + enddo + enddo enddo + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GT/GTpp_self_energy_diag.f90 b/src/GT/GTpp_self_energy_diag.f90 index 6d9e2ed..3612d0e 100644 --- a/src/GT/GTpp_self_energy_diag.f90 +++ b/src/GT/GTpp_self_energy_diag.f90 @@ -1,4 +1,5 @@ -subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) +subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t, & + EcGM,Sig,Z) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -13,13 +14,13 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Om1s(nVVs),Om1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2s(nOOs),Om2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -32,19 +33,31 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho double precision,intent(inout) :: Sig(nBas) double precision,intent(inout) :: Z(nBas) -!---------------------------------------------- -! Occupied part of the T-matrix self-energy -!---------------------------------------------- +! Initialization + + Sig(:) = 0d0 + Z(:) = 0d0 + EcGM = 0d0 + +!--------------------------------------! +! Occupied part of the Tpp self-energy ! +!--------------------------------------! do p=nC+1,nBas-nR do i=nC+1,nO - do cd=1,nVV - eps = e(p) + e(i) - Om1(cd) - num = rho1(p,i,cd)**2 + do cd=1,nVVs + eps = e(p) + e(i) - Om1s(cd) + num = rho1s(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 + do cd=1,nVVt + eps = e(p) + e(i) - Om1t(cd) + num = rho1t(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 @@ -55,14 +68,21 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho do p=nC+1,nBas-nR do a=nO+1,nBas-nR - do kl=1,nOO - eps = e(p) + e(a) - Om2(kl) - num = rho2(p,a,kl)**2 + do kl=1,nOOs + eps = e(p) + e(a) - Om2s(kl) + num = rho2s(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 + + do kl=1,nOOt + eps = e(p) + e(a) - Om2t(kl) + num = rho2t(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 @@ -72,26 +92,40 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho do i=nC+1,nO do j=nC+1,nO - do cd=1,nVV - eps = e(i) + e(j) - Om1(cd) - num = rho1(i,j,cd)**2 + do cd=1,nVVs + eps = e(i) + e(j) - Om1s(cd) + num = rho1s(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + + do cd=1,nVVt + eps = e(i) + e(j) - Om1t(cd) + num = rho1t(i,j,cd)**2 + EcGM = EcGM + num*eps/(eps**2 + eta**2) + enddo + enddo enddo do a=nO+1,nBas-nR do b=nO+1,nBas-nR - do kl=1,nOO - eps = e(a) + e(b) - Om2(kl) - num = rho2(a,b,kl)**2 + do kl=1,nOOs + eps = e(a) + e(b) - Om2s(kl) + num = rho2s(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + + do kl=1,nOOt + eps = e(a) + e(b) - Om2t(kl) + num = rho2t(a,b,kl)**2 + EcGM = EcGM - num*eps/(eps**2 + eta**2) + enddo + enddo enddo + Z(:) = 1d0/(1d0 - Z(:)) + end subroutine diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index 59ed8d8..1f1fdb6 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -1,5 +1,5 @@ -subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & - singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,PHF,cHF,eHF) +subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,PHF,cHF,eHF) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -21,6 +21,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T logical,intent(in) :: dTDA logical,intent(in) :: singlet logical,intent(in) :: triplet + logical,intent(in) :: linearize double precision,intent(in) :: eta logical,intent(in) :: regularize @@ -92,11 +93,11 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & 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),Sig(nBas), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & + eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas), & error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -161,24 +162,15 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- - - EcGM = 0d0 - Sig(:) = 0d0 - Z(:) = 0d0 - - iblock = 3 + iblock = 3 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,Sig,Z) - - iblock = 4 - + iblock = 4 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,Sig,Z) - - Z(:) = 1d0/(1d0 - Z(:)) + call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Solve the quasi-particle equation @@ -186,7 +178,22 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Solve the quasi-particle equation !---------------------------------------------- - eGT(:) = eHF(:) + Sig(:) + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGT(:) = eHF(:) + Sig(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,eOld,eGT) + + end if ! Convergence criteria diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 334307b..0974eb0 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -138,14 +138,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,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),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(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), & - Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & - Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) ! Initialization @@ -217,44 +217,14 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T ! Compute correlation part of the self-energy - EcGM = 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) - if(regularize) then - - 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,Sig,Z) - - end if - iblock = 4 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - if(regularize) then - - 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,Sig,Z) - - end if - - Z(:) = 1d0/(1d0 - Z(:)) + call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis