From b5ccfb734ad4d1042b8da35d44c66fd8826a0546 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 22 Jul 2023 21:40:19 +0200 Subject: [PATCH] remove Vxc and SigX and rename routines in UGT --- src/GT/G0T0eh.f90 | 10 +- src/GT/G0T0pp.f90 | 14 +-- src/GT/{UG0T0.f90 => UG0T0pp.f90} | 97 ++++++------------- ...atrix.f90 => UGTpp_excitation_density.f90} | 4 +- ...x.f90 => UGTpp_renormalization_factor.f90} | 29 +++--- ...ergy_Tmatrix.f90 => UGTpp_self_energy.f90} | 44 ++++----- ...ix_diag.f90 => UGTpp_self_energy_diag.f90} | 44 ++++----- src/GT/evGTeh.f90 | 16 +-- src/GT/evGTpp.f90 | 12 +-- src/GT/{evUGT.f90 => evUGTpp.f90} | 77 +++++---------- src/GT/{qsUGT.f90 => qsUGTpp.f90} | 58 ++++------- src/GW/G0W0.f90 | 12 +-- src/GW/QP_graph.f90 | 6 +- src/GW/UG0W0.f90 | 14 +-- src/GW/evGW.f90 | 14 +-- src/GW/evUGW.f90 | 17 +--- src/GW/unrestricted_QP_graph.f90 | 6 +- src/QuAcK/QuAcK.f90 | 47 +++++---- 18 files changed, 187 insertions(+), 334 deletions(-) rename src/GT/{UG0T0.f90 => UG0T0pp.f90} (61%) rename src/GT/{unrestricted_excitation_density_Tmatrix.f90 => UGTpp_excitation_density.f90} (96%) rename src/GT/{unrestricted_renormalization_factor_Tmatrix.f90 => UGTpp_renormalization_factor.f90} (69%) rename src/GT/{unrestricted_self_energy_Tmatrix.f90 => UGTpp_self_energy.f90} (78%) rename src/GT/{unrestricted_self_energy_Tmatrix_diag.f90 => UGTpp_self_energy_diag.f90} (76%) rename src/GT/{evUGT.f90 => evUGTpp.f90} (72%) rename src/GT/{qsUGT.f90 => qsUGTpp.f90} (82%) diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 8b3a4dc..fc7f1a2 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,6 +1,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) + ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) ! Perform ehG0T0 calculation @@ -37,7 +37,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) @@ -53,7 +52,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, double precision :: EcGM double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) - double precision,allocatable :: SigX(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: Om(:) @@ -99,7 +97,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),SigX(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & + allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),eGT(nBas),eGTlin(nBas)) !-------------------! @@ -123,8 +121,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, ! Compute GW self-energy ! !------------------------! - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - if(regularize) then write(*,*) 'Regularization not yet implemented at the G0T0eh level!' @@ -140,7 +136,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, ! Solve the quasi-particle equation ! !-----------------------------------! - eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) + eGTlin(:) = eHF(:) + Z(:)*Sig(:) ! Linearized or graphical solution? diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index ee590ec..e1c5c6c 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -1,5 +1,5 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -31,7 +31,6 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) @@ -59,7 +58,6 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp double precision,allocatable :: X2ab(:,:),X2aa(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:) - double precision,allocatable :: SigX(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: eGT(:) @@ -91,7 +89,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp 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),Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas)) + Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas)) !---------------------------------------------- ! alpha-beta block @@ -177,17 +175,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp Z(:) = 1d0/(1d0 - Z(:)) -!---------------------------------------------- -! Compute the exchange part of the self-energy -!---------------------------------------------- - - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - !---------------------------------------------- ! Solve the quasi-particle equation !---------------------------------------------- - eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) + eGTlin(:) = eHF(:) + Z(:)*Sig(:) if(linearize) then diff --git a/src/GT/UG0T0.f90 b/src/GT/UG0T0pp.f90 similarity index 61% rename from src/GT/UG0T0.f90 rename to src/GT/UG0T0pp.f90 index 5b1fee5..f6fbb86 100644 --- a/src/GT/UG0T0.f90 +++ b/src/GT/UG0T0pp.f90 @@ -1,7 +1,7 @@ -subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & - spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & - nR,nS,ENuc,EUHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) +subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & + spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & + nR,nS,ENuc,EUHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -32,7 +32,6 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & integer,intent(in) :: nS(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF - double precision,intent(in) :: Vxc(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin) @@ -61,14 +60,12 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) - double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigT(:,:) double precision,allocatable :: Z(:,:) + double precision,allocatable :: eG0T0(:,:) ! Output variables - double precision,intent(out) :: eG0T0(nBas,nspin) - ! Hello world write(*,*) @@ -106,7 +103,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & - SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin)) + SigT(nBas,nspin),Z(nBas,nspin), & + eG0T0(nBas,nspin)) !---------------------------------------------- ! alpha-beta block @@ -118,10 +116,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! 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)) + 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) @@ -137,10 +133,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! Compute linear response - 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)) + 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) @@ -157,10 +151,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! Compute linear response - 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)) + 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) @@ -180,63 +172,42 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & iblock = 3 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & - rho1ab,X2ab,Y2ab,rho2ab) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & + rho1ab,X2ab,Y2ab,rho2ab) !alpha-alpha block iblock = 4 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & - rho1aa,X2aa,Y2aa,rho2aa) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & + rho1aa,X2aa,Y2aa,rho2aa) !beta-beta block iblock = 7 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & - rho1bb,X2bb,Y2bb,rho2bb) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & + rho1bb,X2bb,Y2bb,rho2bb) - call unrestricted_self_energy_Tmatrix_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_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 unrestricted_renormalization_factor_Tmatrix(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) + call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,& + Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z) Z(:,:) = 1d0/(1d0 - Z(:,:)) -!---------------------------------------------- -! Compute the exchange part of the self-energy -!---------------------------------------------- - - do is=1,nspin - call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI,SigX(:,is)) - end do -!call matout(nBas,nspin,SigX) !---------------------------------------------- ! Solve the quasi-particle equation !---------------------------------------------- if(linearize) then -! eG0T0(:) = eHF(:) + Z(:)*SigT(:) - eG0T0(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigT(:,:) - Vxc(:,:)) + eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:) -! call matout(nBas,1,SigX) -! call matout(nBas,1,Vxc) -! call matout(nBas,1,eG0T0(:,1)*HaToeV) -! call matout(nBas,nspin,SigT*HaToeV) else - eG0T0(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) + eG0T0(:,:) = eHF(:,:) + SigT(:,:) end if @@ -251,20 +222,16 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ispin = 1 iblock = 3 - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & - Om2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) !alpha-alpha block ispin = 2 iblock = 4 - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & - Om2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) Ecaa = EcRPA(2) @@ -272,10 +239,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & iblock = 7 - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & - Om2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin)) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb diff --git a/src/GT/unrestricted_excitation_density_Tmatrix.f90 b/src/GT/UGTpp_excitation_density.f90 similarity index 96% rename from src/GT/unrestricted_excitation_density_Tmatrix.f90 rename to src/GT/UGTpp_excitation_density.f90 index 08b0f3c..a4bccb5 100644 --- a/src/GT/unrestricted_excitation_density_Tmatrix.f90 +++ b/src/GT/UGTpp_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nH,nP,ERI_aaaa,ERI_aabb,ERI_bbbb,X1,Y1,rho1,X2,Y2,rho2) +subroutine UGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nH,nP,ERI_aaaa,ERI_aabb,ERI_bbbb,X1,Y1,rho1,X2,Y2,rho2) ! Compute excitation densities for T-matrix self-energy @@ -219,4 +219,4 @@ subroutine unrestricted_excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nH,nP, end if -end subroutine unrestricted_excitation_density_Tmatrix +end subroutine diff --git a/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 b/src/GT/UGTpp_renormalization_factor.f90 similarity index 69% rename from src/GT/unrestricted_renormalization_factor_Tmatrix.f90 rename to src/GT/UGTpp_renormalization_factor.f90 index 918f22d..a66e1fc 100644 --- a/src/GT/unrestricted_renormalization_factor_Tmatrix.f90 +++ b/src/GT/UGTpp_renormalization_factor.f90 @@ -1,8 +1,5 @@ -subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, & - nPaa,nPab,nPbb,e,Omega1aa,Omega1ab, & - Omega1bb,rho1aa,rho1ab,rho1bb, & - Omega2aa,Omega2ab,Omega2bb,rho2aa, & - rho2ab,rho2bb,Z) +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 @@ -16,10 +13,10 @@ subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa integer,intent(in) :: nHaa,nHab,nHbb integer,intent(in) :: nPaa,nPab,nPbb double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega1aa(nPaa),Omega1ab(nPab),Omega1bb(nPbb) + 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) :: Omega2aa(nHaa),Omega2ab(nHab),Omega2bb(nHbb) + 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) @@ -39,14 +36,14 @@ subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa 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) - Omega1aa(cd) + 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) - Omega1ab(cd) + 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 @@ -57,14 +54,14 @@ subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa 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) - Omega2aa(kl) + 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) - Omega2ab(kl) + 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 @@ -77,14 +74,14 @@ subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa 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) - Omega1bb(cd) + 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) - Omega1ab(cd) + 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 @@ -95,17 +92,17 @@ subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa 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) - Omega2bb(kl) + 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) - Omega2ab(kl) + 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 unrestricted_renormalization_factor_Tmatrix +end subroutine diff --git a/src/GT/unrestricted_self_energy_Tmatrix.f90 b/src/GT/UGTpp_self_energy.f90 similarity index 78% rename from src/GT/unrestricted_self_energy_Tmatrix.f90 rename to src/GT/UGTpp_self_energy.f90 index 3f386dd..ba2fd06 100644 --- a/src/GT/unrestricted_self_energy_Tmatrix.f90 +++ b/src/GT/UGTpp_self_energy.f90 @@ -1,7 +1,5 @@ -subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& - nPab,nPbb,e,Omega1aa,Omega1ab,Omega1bb,& - rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& - Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) +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) ! Compute the correlation part of the T-matrix self-energy @@ -19,10 +17,10 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, integer,intent(in) :: nHaa,nHab,nHbb integer,intent(in) :: nPaa,nPab,nPbb double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega1aa(nPaa),Omega1ab(nPab),Omega1bb(nPbb) + 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) :: Omega2aa(nHaa),Omega2ab(nHab),Omega2bb(nHbb) + 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) @@ -46,14 +44,14 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do q=nC(1)+1,nBas-nR(1) do i=nC(1)+1,nO(1) do cd=1,nPaa - eps = e(p,1) + e(i,1) - Omega1aa(cd) + eps = e(p,1) + e(i,1) - Om1aa(cd) SigT(p,q,1) = SigT(p,q,1) + rho1aa(p,i,cd)*rho1aa(q,i,cd)*eps/(eps**2 + eta**2) enddo enddo do i=nC(2)+1,nO(2) do cd=1,nPab - eps = e(p,1) + e(i,1) - Omega1ab(cd) + eps = e(p,1) + e(i,1) - Om1ab(cd) SigT(p,q,1) = SigT(p,q,1) + rho1ab(p,i,cd)*rho1ab(q,i,cd)*eps/(eps**2 + eta**2) end do end do @@ -66,14 +64,14 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do q=nC(2)+1,nBas-nR(2) do i=nC(2)+1,nO(2) do cd=1,nPbb - eps = e(p,2) + e(i,2) - Omega1bb(cd) + eps = e(p,2) + e(i,2) - Om1bb(cd) SigT(p,q,2) = SigT(p,q,2) + rho1bb(p,i,cd)*rho1bb(q,i,cd)*eps/(eps**2 + eta**2) enddo enddo do i=nC(2)+1,nO(2) do cd=1,nPab - eps = e(p,2) + e(i,2) - Omega1ab(cd) + eps = e(p,2) + e(i,2) - Om1ab(cd) SigT(p,q,2) = SigT(p,q,2) + rho1ab(p,i,cd)*rho1ab(q,i,cd)*eps/(eps**2 + eta**2) end do end do @@ -90,14 +88,14 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do q=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) - Omega2aa(kl) + eps = e(p,1) + e(a,1) - Om2aa(kl) SigT(p,q,1) = SigT(p,q,1) + rho2aa(p,a,kl)*rho2aa(q,a,kl)*eps/(eps**2 + eta**2) enddo end do do a=nO(1)+1,nBas-nR(1) do kl=1,nHab - eps = e(p,1) + e(a,1) - Omega2ab(kl) + eps = e(p,1) + e(a,1) - Om2ab(kl) SigT(p,q,1) = SigT(p,q,1) + rho2ab(p,a,kl)*rho2ab(q,a,kl)*eps/(eps**2 + eta**2) end do end do @@ -110,14 +108,14 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do q=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) - Omega2bb(kl) + eps = e(p,2) + e(a,2) - Om2bb(kl) SigT(p,q,2) = SigT(p,q,2) + rho2bb(p,a,kl)*rho2bb(q,a,kl)*eps/(eps**2 + eta**2) enddo enddo do a=nO(2)+1,nBas-nR(2) do kl=1,nHab - eps = e(p,2) + e(a,2) - Omega2ab(kl) + eps = e(p,2) + e(a,2) - Om2ab(kl) SigT(p,q,2) = SigT(p,q,2) + rho2ab(p,a,kl)*rho2ab(q,a,kl)*eps/(eps**2 + eta**2) end do end do @@ -133,7 +131,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, 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) - Omega1aa(cd) + eps = e(i,1) + e(j,1) - Om1aa(cd) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2) enddo enddo @@ -142,7 +140,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do i=nC(1)+1,nO(1) do j=nC(2)+1,nO(2) do cd=1,nPab - eps = e(i,1) + e(j,1) - Omega1ab(cd) + eps = e(i,1) + e(j,1) - Om1ab(cd) EcGM(1) = EcGM(1) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) end do end do @@ -151,7 +149,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do a=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1) do kl=1,nHaa - eps = e(a,1) + e(b,1) - Omega2aa(kl) + eps = e(a,1) + e(b,1) - Om2aa(kl) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -160,7 +158,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, 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) - Omega2ab(kl) + eps = e(a,1) + e(b,1) - Om2ab(kl) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -171,7 +169,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do i=nC(2)+1,nO(2) do j=nC(2)+1,nO(2) do cd=1,nPbb - eps = e(i,2) + e(j,2) - Omega1bb(cd) + eps = e(i,2) + e(j,2) - Om1bb(cd) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2) enddo enddo @@ -180,7 +178,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do i=nC(1)+1,nO(1) do j=nC(2)+1,nO(2) do cd=1,nPab - eps = e(i,2) + e(j,2) - Omega1ab(cd) + eps = e(i,2) + e(j,2) - Om1ab(cd) EcGM(2) = EcGM(2) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) end do end do @@ -189,7 +187,7 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, do a=nO(1)+1,nBas-nR(1) do b=nO(2)+1,nBas-nR(2) do kl=1,nHab - eps = e(a,2) + e(b,2) - Omega2ab(kl) + eps = e(a,2) + e(b,2) - Om2ab(kl) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -198,10 +196,10 @@ subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, 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) - Omega2bb(kl) + eps = e(a,2) + e(b,2) - Om2bb(kl) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo enddo -end subroutine unrestricted_self_energy_Tmatrix +end subroutine diff --git a/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 b/src/GT/UGTpp_self_energy_diag.f90 similarity index 76% rename from src/GT/unrestricted_self_energy_Tmatrix_diag.f90 rename to src/GT/UGTpp_self_energy_diag.f90 index fea06a4..9582b13 100644 --- a/src/GT/unrestricted_self_energy_Tmatrix_diag.f90 +++ b/src/GT/UGTpp_self_energy_diag.f90 @@ -1,7 +1,5 @@ -subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& - nPab,nPbb,e,Omega1aa,Omega1ab,Omega1bb,& - rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& - Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) +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) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -19,10 +17,10 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, integer,intent(in) :: nHaa,nHab,nHbb integer,intent(in) :: nPaa,nPab,nPbb double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega1aa(nPaa),Omega1ab(nPab),Omega1bb(nPbb) + 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) :: Omega2aa(nHaa),Omega2ab(nHab),Omega2bb(nHbb) + 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) @@ -45,14 +43,14 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega1aa(cd) + eps = e(p,1) + e(i,1) - Om1aa(cd) SigT(p,1) = SigT(p,1) + rho1aa(p,i,cd)**2*eps/(eps**2 + eta**2) enddo enddo do i=nC(2)+1,nO(2) do cd=1,nPab - eps = e(p,1) + e(i,1) - Omega1ab(cd) + eps = e(p,1) + e(i,1) - Om1ab(cd) SigT(p,1) = SigT(p,1) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) end do end do @@ -63,14 +61,14 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega1bb(cd) + eps = e(p,2) + e(i,2) - Om1bb(cd) SigT(p,2) = SigT(p,2) + rho1bb(p,i,cd)**2*eps/(eps**2 + eta**2) enddo enddo do i=nC(2)+1,nO(2) do cd=1,nPab - eps = e(p,2) + e(i,2) - Omega1ab(cd) + eps = e(p,2) + e(i,2) - Om1ab(cd) SigT(p,2) = SigT(p,2) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2) end do end do @@ -85,14 +83,14 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega2aa(kl) + eps = e(p,1) + e(a,1) - Om2aa(kl) SigT(p,1) = SigT(p,1) + rho2aa(p,a,kl)**2*eps/(eps**2 + eta**2) enddo end do do a=nO(1)+1,nBas-nR(1) do kl=1,nHab - eps = e(p,1) + e(a,1) - Omega2ab(kl) + eps = e(p,1) + e(a,1) - Om2ab(kl) SigT(p,1) = SigT(p,1) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) end do end do @@ -103,14 +101,14 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega2bb(kl) + eps = e(p,2) + e(a,2) - Om2bb(kl) SigT(p,2) = SigT(p,2) + rho2bb(p,a,kl)**2*eps/(eps**2 + eta**2) enddo enddo do a=nO(2)+1,nBas-nR(2) do kl=1,nHab - eps = e(p,2) + e(a,2) - Omega2ab(kl) + eps = e(p,2) + e(a,2) - Om2ab(kl) SigT(p,2) = SigT(p,2) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2) end do end do @@ -127,7 +125,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega1aa(cd) + eps = e(i,1) + e(j,1) - Om1aa(cd) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2) enddo enddo @@ -136,7 +134,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, do i=nC(1)+1,nO(1) do j=nC(2)+1,nO(2) do cd=1,nPab - eps = e(i,1) + e(j,1) - Omega1ab(cd) + eps = e(i,1) + e(j,1) - Om1ab(cd) EcGM(1) = EcGM(1) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) end do end do @@ -145,7 +143,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, do a=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1) do kl=1,nHaa - eps = e(a,1) + e(b,1) - Omega2aa(kl) + eps = e(a,1) + e(b,1) - Om2aa(kl) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -154,7 +152,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega2ab(kl) + eps = e(a,1) + e(b,1) - Om2ab(kl) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -165,7 +163,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, do i=nC(2)+1,nO(2) do j=nC(2)+1,nO(2) do cd=1,nPbb - eps = e(i,2) + e(j,2) - Omega1bb(cd) + eps = e(i,2) + e(j,2) - Om1bb(cd) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2) enddo enddo @@ -174,7 +172,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, do i=nC(1)+1,nO(1) do j=nC(2)+1,nO(2) do cd=1,nPab - eps = e(i,2) + e(j,2) - Omega1ab(cd) + eps = e(i,2) + e(j,2) - Om1ab(cd) EcGM(2) = EcGM(2) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2) end do end do @@ -183,7 +181,7 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, do a=nO(1)+1,nBas-nR(1) do b=nO(2)+1,nBas-nR(2) do kl=1,nHab - eps = e(a,2) + e(b,2) - Omega2ab(kl) + eps = e(a,2) + e(b,2) - Om2ab(kl) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo @@ -192,10 +190,10 @@ subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab, 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) - Omega2bb(kl) + eps = e(a,2) + e(b,2) - Om2bb(kl) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2) enddo enddo enddo -end subroutine unrestricted_self_energy_Tmatrix_diag +end subroutine diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index 5072bb0..049e861 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -1,6 +1,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & - cHF,eHF,Vxc) + cHF,eHF) ! Perform self-consistent eigenvalue-only ehGT calculation @@ -39,7 +39,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -64,7 +63,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d double precision,allocatable :: eGT(:) double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) - double precision,allocatable :: SigX(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) @@ -103,13 +101,9 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),Sig(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & + allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas)) -! Compute the exchange part of the self-energy - - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - ! Initialization nSCF = 0 @@ -155,7 +149,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Solve the quasi-particle equation - eGTlin(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) + eGTlin(:) = eHF(:) + Sig(:) ! Linearized or graphical solution? @@ -170,9 +164,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' ! write(*,*) -! -! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGTlin,eGT,regularize) - + end if ! Convergence criteria diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index acbd3d2..801ed93 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_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) + singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -35,7 +35,6 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -69,7 +68,6 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T double precision,allocatable :: X2s(:,:),X2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) - double precision,allocatable :: SigX(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) @@ -99,13 +97,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T 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),Sig(nBas), & + eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas), & error_diis(nBas,max_diis),e_diis(nBas,max_diis)) -! Compute the exchange part of the self-energy - - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - ! Initialization nSCF = 0 @@ -193,7 +187,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Solve the quasi-particle equation !---------------------------------------------- - eGT(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) + eGT(:) = eHF(:) + Sig(:) ! Convergence criteria diff --git a/src/GT/evUGT.f90 b/src/GT/evUGTpp.f90 similarity index 72% rename from src/GT/evUGT.f90 rename to src/GT/evUGTpp.f90 index bea6bcf..6b20a5d 100644 --- a/src/GT/evUGT.f90 +++ b/src/GT/evUGTpp.f90 @@ -1,8 +1,7 @@ -subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & - TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& - eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO,ERI_aaaa, & - ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF, & - Vxc,eG0T0) +subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & + TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& + eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO,ERI_aaaa, & + ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -34,7 +33,6 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & integer,intent(in) :: nS(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF - double precision,intent(in) :: Vxc(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin) @@ -44,7 +42,6 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) - double precision,intent(in) :: eG0T0(nBas,nspin) ! Local variables integer :: nSCF @@ -67,7 +64,6 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) - double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigT(:,:) double precision,allocatable :: Z(:,:) double precision,allocatable :: eGT(:,:) @@ -113,18 +109,13 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & - SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin), & + SigT(nBas,nspin),Z(nBas,nspin), & eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), & e_diis(nBas,max_diis,nspin)) !---------------------------------------------- ! Compute the exchange part of the self-energy !---------------------------------------------- - - do is=1,nspin - call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI_AO, & - SigX(:,is)) - end do !Initialization @@ -133,7 +124,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & Conv = 1d0 e_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0 - eGT(:,:) = eG0T0(:,:) + eGT(:,:) = eHF(:,:) eOld(:,:) = eGT(:,:) Z(:,:) = 1d0 rcond(:) = 0d0 @@ -154,13 +145,9 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & - Om2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) -! EcRPA(ispin) = 1d0*EcRPA(ispin) - !---------------------------------------------- ! alpha-alpha block !---------------------------------------------- @@ -170,14 +157,9 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & - Om2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) -! EcRPA(ispin) = 2d0*EcRPA(ispin) -! EcRPA(ispin) = 3d0*EcRPA(ispin) - !---------------------------------------------- ! beta-beta block !---------------------------------------------- @@ -187,13 +169,8 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,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 ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin)) !---------------------------------------------- ! Compute T-matrix version of the self-energy @@ -207,35 +184,27 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & iblock = 3 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & - rho1ab,X2ab,Y2ab,rho2ab) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & + rho1ab,X2ab,Y2ab,rho2ab) !alpha-alpha block iblock = 4 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & - rho1aa,X2aa,Y2aa,rho2aa) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & + rho1aa,X2aa,Y2aa,rho2aa) !beta-beta block iblock = 7 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & - rho1bb,X2bb,Y2bb,rho2bb) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & + rho1bb,X2bb,Y2bb,rho2bb) - call unrestricted_self_energy_Tmatrix_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_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 unrestricted_renormalization_factor_Tmatrix(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) + 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(:,:)) @@ -244,7 +213,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Solve the quasi-particle equation !---------------------------------------------- - eGT(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) + eGT(:,:) = eHF(:,:) + SigT(:,:) ! Convergence criteria diff --git a/src/GT/qsUGT.f90 b/src/GT/qsUGTpp.f90 similarity index 82% rename from src/GT/qsUGT.f90 rename to src/GT/qsUGTpp.f90 index 9508374..f69054b 100644 --- a/src/GT/qsUGT.f90 +++ b/src/GT/qsUGTpp.f90 @@ -1,7 +1,7 @@ -subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & - TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& - eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,& - ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) +subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & + TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& + eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,& + ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -208,10 +208,8 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & - Om2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) ! EcRPA(ispin) = 1d0*EcRPA(ispin) @@ -224,13 +222,9 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & - Om2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) -! EcRPA(ispin) = 2d0*EcRPA(ispin) -! EcRPA(ispin) = 3d0*EcRPA(ispin) Ecaa = EcRPA(2) !---------------------------------------------- @@ -242,13 +236,9 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & - Om2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin)) -! EcRPA(ispin) = 2d0*EcRPA(ispin) -! EcRPA(ispin) = 3d0*EcRPA(ispin) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb EcRPA(1) = EcRPA(1) - EcRPA(2) @@ -266,35 +256,27 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & iblock = 3 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & - rho1ab,X2ab,Y2ab,rho2ab) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & + rho1ab,X2ab,Y2ab,rho2ab) !alpha-alpha block iblock = 4 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & - rho1aa,X2aa,Y2aa,rho2aa) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & + rho1aa,X2aa,Y2aa,rho2aa) !beta-beta block iblock = 7 - call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & - ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & - rho1bb,X2bb,Y2bb,rho2bb) + call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & + rho1bb,X2bb,Y2bb,rho2bb) - call unrestricted_self_energy_Tmatrix(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_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 unrestricted_renormalization_factor_Tmatrix(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) + 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(:,:)) diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 0c86e37..5ae8a1b 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,5 +1,5 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) ! Perform G0W0 calculation @@ -36,7 +36,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) @@ -51,7 +50,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT double precision :: EcGM double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) - double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) double precision,allocatable :: Om(:) @@ -97,7 +95,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),SigX(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & + allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & eGW(nBas),eGWlin(nBas)) !-------------------! @@ -121,8 +119,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT ! Compute GW self-energy ! !------------------------! - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - if(regularize) then call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC) @@ -138,7 +134,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT ! Solve the quasi-particle equation ! !-----------------------------------! - eGWlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) + eGWlin(:) = eHF(:) + Z(:)*SigC(:) ! Linearized or graphical solution? @@ -154,7 +150,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGWlin,eGW,regularize) + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize) end if diff --git a/src/GW/QP_graph.f90 b/src/GW/QP_graph.f90 index 93faf9d..7c5bed2 100644 --- a/src/GW/QP_graph.f90 +++ b/src/GW/QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW,regularize) +subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize) ! Compute the graphical solution of the QP equation @@ -15,8 +15,6 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW,re integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: SigX(nBas) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) @@ -57,7 +55,7 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW,re sigC = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize) dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize) - f = w - eHF(p) - SigX(p) - sigC + Vxc(p) + f = w - eHF(p) - SigC df = 1d0 - dsigC w = w - f/df diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index f3994f6..20aa140 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -1,6 +1,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW) + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform unrestricted G0W0 calculation @@ -42,7 +42,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin) - double precision,intent(in) :: Vxc(nBas,nspin) ! Local variables @@ -53,7 +52,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons double precision :: EcGM(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) - double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) integer :: nS_aa,nS_bb,nS_sc @@ -101,7 +99,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(SigX(nBas,nspin),SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin), & + allocate(SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin), & OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) !-------------------! @@ -127,10 +125,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Compute self-energy and renormalization factor ! !------------------------------------------------! - do is=1,nspin - call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI,SigX(:,is)) - end do - if(regularize) then call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) @@ -147,7 +141,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Solve the quasi-particle equation ! !-----------------------------------! - eGWlin(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigC(:,:) - Vxc(:,:)) + eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) if(linearize) then @@ -161,7 +155,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Find graphical solution of the QP equation do is=1,nspin - call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),SigX(:,is),Vxc(:,is), & + call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), & OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) end do diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 5fd8178..9f17631 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,6 +1,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & - cHF,eHF,Vxc) + cHF,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -39,7 +39,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -64,7 +63,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop double precision,allocatable :: eGW(:) double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) - double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) @@ -100,13 +98,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas), & + allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), & Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) -! Compute the exchange part of the self-energy - - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - ! Initialization nSCF = 0 @@ -152,7 +146,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Solve the quasi-particle equation - eGW(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + eGW(:) = eHF(:) + SigC(:) ! Linearized or graphical solution? @@ -168,7 +162,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGW,eGW,regularize) + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize) end if diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 248926a..19690e7 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & dBSE,dTDA,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & - EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) + EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -37,8 +37,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) - double precision,intent(in) :: Vxc(nBas,nspin) - double precision,intent(in) :: eG0W0(nBas,nspin) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) @@ -67,7 +65,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,allocatable :: eOld(:,:) double precision,allocatable :: Z(:,:) integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: SigX(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) @@ -107,16 +104,10 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigX(nBas,nspin),SigC(nBas,nspin), & + allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin), & OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), & error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) -! Compute the exchange part of the self-energy - - do is=1,nspin - call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI_AO,SigX(:,is)) - end do - ! Initialization nSCF = 0 @@ -125,7 +116,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, Conv = 1d0 e_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0 - eGW(:,:) = eG0W0(:,:) + eGW(:,:) = eHF(:,:) eOld(:,:) = eGW(:,:) Z(:,:) = 1d0 rcond(:) = 0d0 @@ -167,7 +158,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Solve the quasi-particle equation ! !-----------------------------------! - eGW(:,:) = eHF(:,:) + SigX(:,:) + SigC(:,:) - Vxc(:,:) + eGW(:,:) = eHF(:,:) + SigC(:,:) ! Convergence criteria diff --git a/src/GW/unrestricted_QP_graph.f90 b/src/GW/unrestricted_QP_graph.f90 index 14b7653..6162b7c 100644 --- a/src/GW/unrestricted_QP_graph.f90 +++ b/src/GW/unrestricted_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW) +subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) ! Compute the graphical solution of the QP equation @@ -15,8 +15,6 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho, integer,intent(in) :: nS double precision,intent(in) :: eta double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: SigX(nBas) - double precision,intent(in) :: Vxc(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS,nspin) @@ -56,7 +54,7 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho, sigC = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) dsigC = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) - f = w - eHF(p) - SigX(p) + Vxc(p) - sigC + f = w - eHF(p) - sigC df = 1d0 - dsigC w = w - f/df diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index a7fcb25..dbd8600 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -28,7 +28,6 @@ program QuAcK double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:) - double precision,allocatable :: Vxc(:,:) logical :: doACFDT logical :: exchange_kernel @@ -186,7 +185,7 @@ program QuAcK allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), & V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), & - dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin),F_AO(nBas,nBas)) + dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas)) ! Read integrals @@ -217,7 +216,7 @@ program QuAcK call wall_time(start_HF) call HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,n_diis_HF, & guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, & - ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF,Vxc) + ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -808,11 +807,11 @@ program QuAcK call UG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF,Vxc) + dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF) end if call cpu_time(end_GW) @@ -835,13 +834,13 @@ program QuAcK call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA, & dBSE,dTDA,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc, & EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & - PHF,cHF,epsHF,Vxc) + PHF,cHF,epsHF) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,linGW,eta_GW,regGW, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF) end if call cpu_time(end_GW) @@ -952,15 +951,15 @@ program QuAcK if(unrestricted) then - call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA, & - spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & - nR,nS,ENuc,EHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & - dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF,Vxc) + call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA, & + spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & + nR,nS,ENuc,EHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) + linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF) end if @@ -982,18 +981,18 @@ program QuAcK if(unrestricted) then - call evUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& - eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, & - dipole_int_bb,PHF,cHF,epsHF,Vxc) + call evUGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& + eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, & + dipole_int_bb,PHF,cHF,epsHF) else call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO, & - PHF,cHF,epsHF,Vxc) + PHF,cHF,epsHF) end if @@ -1015,10 +1014,10 @@ program QuAcK if(unrestricted) then - call qsUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T, & - TDA,dBSE,dTDA,spin_conserved,spin_flip,eta_GT,regGT,nBas,nC,nO,nV,& - nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,& - ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + call qsUGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T, & + TDA,dBSE,dTDA,spin_conserved,spin_flip,eta_GT,regGT,nBas,nC,nO,nV,& + nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,& + ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else call qsGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & @@ -1051,7 +1050,7 @@ program QuAcK else call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) + linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF) end if @@ -1076,7 +1075,7 @@ program QuAcK call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet,linGT,eta_GT,regGT, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF) end if call cpu_time(end_GT)