LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-09 12:44:09 +01:00
Code Releases Activity

remove Vxc and SigX and rename routines in UGT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-22 21:40:19 +02:00
parent 1ea1d49a4e
commit b5ccfb734a
18 changed files with 187 additions and 334 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/GT/G0T0eh.f90
View File

@ -1,6 +1,6 @@
subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 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, & 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 ! 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_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(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) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(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 :: EcGM
double precision,allocatable :: Aph(:,:) double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:) double precision,allocatable :: Bph(:,:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: Sig(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: Om(:) double precision,allocatable :: Om(:)
@ -99,7 +97,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
! Memory allocation ! 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)) 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 ! ! Compute GW self-energy !
!------------------------! !------------------------!
call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
if(regularize) then if(regularize) then
write(*,*) 'Regularization not yet implemented at the G0T0eh level!' 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 ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) eGTlin(:) = eHF(:) + Z(:)*Sig(:)
! Linearized or graphical solution? ! Linearized or graphical solution?

14
src/GT/G0T0pp.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 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) ! 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 integer,intent(in) :: nS
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(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 :: X2ab(:,:),X2aa(:,:)
double precision,allocatable :: Y2ab(:,:),Y2aa(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:)
double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: Sig(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: eGT(:) 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), & Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), &
Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), &
rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,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 ! alpha-beta block
@ -177,17 +175,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
Z(:) = 1d0/(1d0 - Z(:)) 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 ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + Sig(:) - Vxc(:)) eGTlin(:) = eHF(:) + Z(:)*Sig(:)
if(linearize) then if(linearize) then

97
src/GT/UG0T0.f90 → src/GT/UG0T0pp.f90
View File

@ -1,7 +1,7 @@
subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, &
nR,nS,ENuc,EUHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & nR,nS,ENuc,EUHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! 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) integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
double precision,intent(in) :: Vxc(nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(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 :: X2ab(:,:),X2aa(:,:),X2bb(:,:)
double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:)
double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:)
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigT(:,:) double precision,allocatable :: SigT(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
double precision,allocatable :: eG0T0(:,:)
! Output variables ! Output variables
double precision,intent(out) :: eG0T0(nBas,nspin)
! Hello world ! Hello world
write(*,*) 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), & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), &
Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), &
rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin)) SigT(nBas,nspin),Z(nBas,nspin), &
eG0T0(nBas,nspin))
!---------------------------------------------- !----------------------------------------------
! alpha-beta block ! alpha-beta block
@ -118,10 +116,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, &
nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, &
Om2ab,X2ab,Y2ab,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin) ! EcRPA(ispin) = 1d0*EcRPA(ispin)
@ -137,10 +133,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, &
nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, &
Om2aa,X2aa,Y2aa,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*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 ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, &
nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, &
Om2bb,X2bb,Y2bb,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*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 iblock = 3
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & rho1ab,X2ab,Y2ab,rho2ab)
rho1ab,X2ab,Y2ab,rho2ab)
!alpha-alpha block !alpha-alpha block
iblock = 4 iblock = 4
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & rho1aa,X2aa,Y2aa,rho2aa)
rho1aa,X2aa,Y2aa,rho2aa)
!beta-beta block !beta-beta block
iblock = 7 iblock = 7
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & rho1bb,X2bb,Y2bb,rho2bb)
rho1bb,X2bb,Y2bb,rho2bb)
call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,Om1bb,&
nPab,nPbb,eHF,Om1aa,Om1ab,Om1bb,& rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,&
Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,&
nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,& Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Om1bb,rho1aa,rho1ab,rho1bb, &
Om2aa,Om2ab,Om2bb,rho2aa, &
rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - 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 ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
if(linearize) then if(linearize) then
! eG0T0(:) = eHF(:) + Z(:)*SigT(:) eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
eG0T0(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigT(:,:) - Vxc(:,:))
! call matout(nBas,1,SigX)
! call matout(nBas,1,Vxc)
! call matout(nBas,1,eG0T0(:,1)*HaToeV)
! call matout(nBas,nspin,SigT*HaToeV)
else else
eG0T0(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) eG0T0(:,:) = eHF(:,:) + SigT(:,:)
end if end if
@ -251,20 +222,16 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
ispin = 1 ispin = 1
iblock = 3 iblock = 3
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, &
nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, &
Om2ab,X2ab,Y2ab,EcRPA(ispin))
!alpha-alpha block !alpha-alpha block
ispin = 2 ispin = 2
iblock = 4 iblock = 4
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, &
nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, &
Om2aa,X2aa,Y2aa,EcRPA(ispin))
Ecaa = EcRPA(2) Ecaa = EcRPA(2)
@ -272,10 +239,8 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
iblock = 7 iblock = 7
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, &
nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, &
Om2bb,X2bb,Y2bb,EcRPA(ispin))
Ecbb = EcRPA(2) Ecbb = EcRPA(2)
EcRPA(2) = Ecaa + Ecbb EcRPA(2) = Ecaa + Ecbb

4
src/GT/unrestricted_excitation_density_Tmatrix.f90 → src/GT/UGTpp_excitation_density.f90
View File

@ -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 ! 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 if
end subroutine unrestricted_excitation_density_Tmatrix end subroutine

29
src/GT/unrestricted_renormalization_factor_Tmatrix.f90 → src/GT/UGTpp_renormalization_factor.f90
View File

@ -1,8 +1,5 @@
subroutine unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb, & subroutine UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab, &
nPaa,nPab,nPbb,e,Omega1aa,Omega1ab, & Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Omega1bb,rho1aa,rho1ab,rho1bb, &
Omega2aa,Omega2ab,Omega2bb,rho2aa, &
rho2ab,rho2bb,Z)
! Compute renormalization factor of the T-matrix self-energy ! 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) :: nHaa,nHab,nHbb
integer,intent(in) :: nPaa,nPab,nPbb integer,intent(in) :: nPaa,nPab,nPbb
double precision,intent(in) :: e(nBas,nspin) 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) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab)
double precision,intent(in) :: rho1bb(nBas,nBas,nPbb) 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) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab)
double precision,intent(in) :: rho2bb(nBas,nBas,nHbb) 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 p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(p,1) + e(i,1) - 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 Z(p,1) = Z(p,1) - rho1aa(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
enddo enddo
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPab 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 Z(p,1) = Z(p,1) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
end do end do
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 p=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa do kl=1,nHaa
eps = e(p,1) + e(a,1) - 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 Z(p,1) = Z(p,1) - rho2aa(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
enddo enddo
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab do kl=1,nHab
eps = e(p,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 Z(p,1) = Z(p,1) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
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 p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPbb do cd=1,nPbb
eps = e(p,2) + e(i,2) - 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 Z(p,2) = Z(p,2) - rho1bb(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
enddo enddo
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPab 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 Z(p,2) = Z(p,2) - rho1ab(p,i,cd)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
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 p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(p,2) + e(a,2) - 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 Z(p,2) = Z(p,2) - rho2bb(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
enddo enddo
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(p,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 Z(p,2) = Z(p,2) - rho2ab(p,a,kl)**2*(eps/(eps**2 + eta**2))**2
enddo enddo
enddo enddo
enddo enddo
end subroutine unrestricted_renormalization_factor_Tmatrix end subroutine

44
src/GT/unrestricted_self_energy_Tmatrix.f90 → src/GT/UGTpp_self_energy.f90
View File

@ -1,7 +1,5 @@
subroutine unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& subroutine UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,&
nPab,nPbb,e,Omega1aa,Omega1ab,Omega1bb,& rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,&
Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
! Compute the correlation part of the T-matrix self-energy ! 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) :: nHaa,nHab,nHbb
integer,intent(in) :: nPaa,nPab,nPbb integer,intent(in) :: nPaa,nPab,nPbb
double precision,intent(in) :: e(nBas,nspin) 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) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab)
double precision,intent(in) :: rho1bb(nBas,nBas,nPbb) 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) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab)
double precision,intent(in) :: rho2bb(nBas,nBas,nHbb) 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 q=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(p,1) + e(i,1) - 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) SigT(p,q,1) = SigT(p,q,1) + rho1aa(p,i,cd)*rho1aa(q,i,cd)*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,1) + e(i,1) - 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) 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
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 q=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPbb do cd=1,nPbb
eps = e(p,2) + e(i,2) - 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) SigT(p,q,2) = SigT(p,q,2) + rho1bb(p,i,cd)*rho1bb(q,i,cd)*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,2) + e(i,2) - 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) 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
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 q=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa do kl=1,nHaa
eps = e(p,1) + e(a,1) - 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) SigT(p,q,1) = SigT(p,q,1) + rho2aa(p,a,kl)*rho2aa(q,a,kl)*eps/(eps**2 + eta**2)
enddo enddo
end do end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(p,1) + e(a,1) - 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) 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
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 q=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(p,2) + e(a,2) - 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) SigT(p,q,2) = SigT(p,q,2) + rho2bb(p,a,kl)*rho2bb(q,a,kl)*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab do kl=1,nHab
eps = e(p,2) + e(a,2) - 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) 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
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 i=nC(1)+1,nO(1)
do j=nC(1)+1,nO(1) do j=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(i,1) + e(j,1) - 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) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPab 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) EcGM(1) = EcGM(1) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2)
end do end do
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 a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHaa 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) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(a,1) + e(b,1) - 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) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(2)+1,nO(2)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPbb 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) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPab 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) EcGM(2) = EcGM(2) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2)
end do end do
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 a=nO(1)+1,nBas-nR(1)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHab 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) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(a,2) + e(b,2) - 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) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo
end subroutine unrestricted_self_energy_Tmatrix end subroutine

44
src/GT/unrestricted_self_energy_Tmatrix_diag.f90 → src/GT/UGTpp_self_energy_diag.f90
View File

@ -1,7 +1,5 @@
subroutine unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& subroutine UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,e,Om1aa,Om1ab,Om1bb,&
nPab,nPbb,e,Omega1aa,Omega1ab,Omega1bb,& rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,&
Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
! Compute diagonal of the correlation part of the T-matrix self-energy ! 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) :: nHaa,nHab,nHbb
integer,intent(in) :: nPaa,nPab,nPbb integer,intent(in) :: nPaa,nPab,nPbb
double precision,intent(in) :: e(nBas,nspin) 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) :: rho1aa(nBas,nBas,nPaa),rho1ab(nBas,nBas,nPab)
double precision,intent(in) :: rho1bb(nBas,nBas,nPbb) 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) :: rho2aa(nBas,nBas,nHaa),rho2ab(nBas,nBas,nHab)
double precision,intent(in) :: rho2bb(nBas,nBas,nHbb) 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 p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(p,1) + e(i,1) - 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) SigT(p,1) = SigT(p,1) + rho1aa(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,1) + e(i,1) - 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) SigT(p,1) = SigT(p,1) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2)
end do end do
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 p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPbb do cd=1,nPbb
eps = e(p,2) + e(i,2) - 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) SigT(p,2) = SigT(p,2) + rho1bb(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do i=nC(2)+1,nO(2) do i=nC(2)+1,nO(2)
do cd=1,nPab do cd=1,nPab
eps = e(p,2) + e(i,2) - 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) SigT(p,2) = SigT(p,2) + rho1ab(p,i,cd)**2*eps/(eps**2 + eta**2)
end do end do
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 p=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHaa do kl=1,nHaa
eps = e(p,1) + e(a,1) - 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) SigT(p,1) = SigT(p,1) + rho2aa(p,a,kl)**2*eps/(eps**2 + eta**2)
enddo enddo
end do end do
do a=nO(1)+1,nBas-nR(1) do a=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(p,1) + e(a,1) - 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) SigT(p,1) = SigT(p,1) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2)
end do end do
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 p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(p,2) + e(a,2) - 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) SigT(p,2) = SigT(p,2) + rho2bb(p,a,kl)**2*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
do kl=1,nHab do kl=1,nHab
eps = e(p,2) + e(a,2) - 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) SigT(p,2) = SigT(p,2) + rho2ab(p,a,kl)**2*eps/(eps**2 + eta**2)
end do end do
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 i=nC(1)+1,nO(1)
do j=nC(1)+1,nO(1) do j=nC(1)+1,nO(1)
do cd=1,nPaa do cd=1,nPaa
eps = e(i,1) + e(j,1) - 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) EcGM(1) = EcGM(1) + rho1aa(i,j,cd)*rho1aa(i,j,cd)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPab 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) EcGM(1) = EcGM(1) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2)
end do end do
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 a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHaa 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) EcGM(1) = EcGM(1) - rho2aa(a,b,kl)*rho2aa(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1) do b=nO(1)+1,nBas-nR(1)
do kl=1,nHab do kl=1,nHab
eps = e(a,1) + e(b,1) - 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) EcGM(1) = EcGM(1) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(2)+1,nO(2)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPbb 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) EcGM(2) = EcGM(2) + rho1bb(i,j,cd)*rho1bb(i,j,cd)*eps/(eps**2 + eta**2)
enddo enddo
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 i=nC(1)+1,nO(1)
do j=nC(2)+1,nO(2) do j=nC(2)+1,nO(2)
do cd=1,nPab 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) EcGM(2) = EcGM(2) + rho1ab(i,j,cd)*rho1ab(i,j,cd)*eps/(eps**2 + eta**2)
end do end do
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 a=nO(1)+1,nBas-nR(1)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHab 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) EcGM(2) = EcGM(2) - rho2ab(a,b,kl)*rho2ab(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
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 a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2) do b=nO(2)+1,nBas-nR(2)
do kl=1,nHbb do kl=1,nHbb
eps = e(a,2) + e(b,2) - 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) EcGM(2) = EcGM(2) - rho2bb(a,b,kl)*rho2bb(a,b,kl)*eps/(eps**2 + eta**2)
enddo enddo
enddo enddo
enddo enddo
end subroutine unrestricted_self_energy_Tmatrix_diag end subroutine

14
src/GT/evGTeh.f90
View File

@ -1,6 +1,6 @@
subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & 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, & 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 ! 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) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,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_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(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) :: 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 :: eGT(:)
double precision,allocatable :: eOld(:) double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: Sig(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Om(:) double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:) double precision,allocatable :: XpY(:,:)
@ -103,13 +101,9 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Memory allocation ! 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)) 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 ! Initialization
nSCF = 0 nSCF = 0
@ -155,7 +149,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Solve the quasi-particle equation ! Solve the quasi-particle equation
eGTlin(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) eGTlin(:) = eHF(:) + Sig(:)
! Linearized or graphical solution? ! Linearized or graphical solution?
@ -170,8 +164,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' ! write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
! write(*,*) ! write(*,*)
!
! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGTlin,eGT,regularize)
end if end if

12
src/GT/evGTpp.f90
View File

@ -1,5 +1,5 @@
subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & 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) ! 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) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,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_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(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) :: 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 :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: Sig(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) 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), & Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,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)) 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 ! Initialization
nSCF = 0 nSCF = 0
@ -193,7 +187,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Solve the quasi-particle equation ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
eGT(:) = eHF(:) + SigX(:) + Sig(:) - Vxc(:) eGT(:) = eHF(:) + Sig(:)
! Convergence criteria ! Convergence criteria

77
src/GT/evUGT.f90 → src/GT/evUGTpp.f90
View File

@ -1,8 +1,7 @@
subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO,ERI_aaaa, & 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, & ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
Vxc,eG0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! 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) integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
double precision,intent(in) :: Vxc(nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(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) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: eG0T0(nBas,nspin)
! Local variables ! Local variables
integer :: nSCF 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 :: X2ab(:,:),X2aa(:,:),X2bb(:,:)
double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:)
double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:)
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigT(:,:) double precision,allocatable :: SigT(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
double precision,allocatable :: eGT(:,:) double precision,allocatable :: eGT(:,:)
@ -113,7 +109,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), &
Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), &
rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
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), & eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), &
e_diis(nBas,max_diis,nspin)) e_diis(nBas,max_diis,nspin))
@ -121,11 +117,6 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute the exchange part of the self-energy ! 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 !Initialization
nSCF = 0 nSCF = 0
@ -133,7 +124,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
Conv = 1d0 Conv = 1d0
e_diis(:,:,:) = 0d0 e_diis(:,:,:) = 0d0
error_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0
eGT(:,:) = eG0T0(:,:) eGT(:,:) = eHF(:,:)
eOld(:,:) = eGT(:,:) eOld(:,:) = eGT(:,:)
Z(:,:) = 1d0 Z(:,:) = 1d0
rcond(:) = 0d0 rcond(:) = 0d0
@ -154,12 +145,8 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, &
nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, &
Om2ab,X2ab,Y2ab,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin)
!---------------------------------------------- !----------------------------------------------
! alpha-alpha block ! alpha-alpha block
@ -170,13 +157,8 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, &
nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, &
Om2aa,X2aa,Y2aa,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
!---------------------------------------------- !----------------------------------------------
! beta-beta block ! beta-beta block
@ -187,13 +169,8 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, &
nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, &
Om2bb,X2bb,Y2bb,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
!---------------------------------------------- !----------------------------------------------
! Compute T-matrix version of the self-energy ! 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 iblock = 3
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & rho1ab,X2ab,Y2ab,rho2ab)
rho1ab,X2ab,Y2ab,rho2ab)
!alpha-alpha block !alpha-alpha block
iblock = 4 iblock = 4
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & rho1aa,X2aa,Y2aa,rho2aa)
rho1aa,X2aa,Y2aa,rho2aa)
!beta-beta block !beta-beta block
iblock = 7 iblock = 7
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & rho1bb,X2bb,Y2bb,rho2bb)
rho1bb,X2bb,Y2bb,rho2bb)
call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& call UGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,&
nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,&
Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,&
nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,& Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Om1bb,rho1aa,rho1ab,rho1bb, &
Om2aa,Om2ab,Om2bb,rho2aa, &
rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - 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 ! Solve the quasi-particle equation
!---------------------------------------------- !----------------------------------------------
eGT(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) eGT(:,:) = eHF(:,:) + SigT(:,:)
! Convergence criteria ! Convergence criteria

58
src/GT/qsUGT.f90 → src/GT/qsUGTpp.f90
View File

@ -1,7 +1,7 @@
subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& 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,& 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) ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GT calculation ! 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 ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, &
nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, &
Om2ab,X2ab,Y2ab,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin) ! EcRPA(ispin) = 1d0*EcRPA(ispin)
@ -224,13 +222,9 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, &
nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, &
Om2aa,X2aa,Y2aa,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
Ecaa = EcRPA(2) Ecaa = EcRPA(2)
!---------------------------------------------- !----------------------------------------------
@ -242,13 +236,9 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Compute linear response ! Compute linear response
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, &
nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, &
Om2bb,X2bb,Y2bb,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin)
Ecbb = EcRPA(2) Ecbb = EcRPA(2)
EcRPA(2) = Ecaa + Ecbb EcRPA(2) = Ecaa + Ecbb
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = EcRPA(1) - EcRPA(2)
@ -266,35 +256,27 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
iblock = 3 iblock = 3
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHab,nPab,ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & rho1ab,X2ab,Y2ab,rho2ab)
rho1ab,X2ab,Y2ab,rho2ab)
!alpha-alpha block !alpha-alpha block
iblock = 4 iblock = 4
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa,ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & rho1aa,X2aa,Y2aa,rho2aa)
rho1aa,X2aa,Y2aa,rho2aa)
!beta-beta block !beta-beta block
iblock = 7 iblock = 7
call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & call UGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb,ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, &
ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & rho1bb,X2bb,Y2bb,rho2bb)
rho1bb,X2bb,Y2bb,rho2bb)
call unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& call UGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,&
nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,&
Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT)
call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& call UGTpp_renormalization_factor(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,&
nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,& Om1bb,rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,Om2bb,rho2aa,rho2ab,rho2bb,Z)
Om1bb,rho1aa,rho1ab,rho1bb, &
Om2aa,Om2ab,Om2bb,rho2aa, &
rho2ab,rho2bb,Z)
Z(:,:) = 1d0/(1d0 - Z(:,:)) Z(:,:) = 1d0/(1d0 - Z(:,:))

12
src/GW/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 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 ! 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_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(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) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(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 :: EcGM
double precision,allocatable :: Aph(:,:) double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:) double precision,allocatable :: Bph(:,:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: Om(:) double precision,allocatable :: Om(:)
@ -97,7 +95,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
! Memory allocation ! 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)) 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 ! ! Compute GW self-energy !
!------------------------! !------------------------!
call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
if(regularize) then if(regularize) then
call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC) 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 ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGWlin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) eGWlin(:) = eHF(:) + Z(:)*SigC(:)
! Linearized or graphical solution? ! 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(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*) 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 end if

6
src/GW/QP_graph.f90
View File

@ -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 ! 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 integer,intent(in) :: nS
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas) 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) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,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) 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) 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 df = 1d0 - dsigC
w = w - f/df w = w - f/df

14
src/GW/UG0W0.f90
View File

@ -1,6 +1,6 @@
subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & 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, & 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 ! 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) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: Vxc(nBas,nspin)
! Local variables ! 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 :: EcGM(nspin)
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin) double precision :: EcAC(nspin)
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigC(:,:) double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc 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_bb = nS(2)
nS_sc = nS_aa + nS_bb 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)) 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 ! ! 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 if(regularize) then
call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) 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 ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGWlin(:,:) = eHF(:,:) + Z(:,:)*(SigX(:,:) + SigC(:,:) - Vxc(:,:)) eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
if(linearize) then 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 ! Find graphical solution of the QP equation
do is=1,nspin 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)) OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do end do

14
src/GW/evGW.f90
View File

@ -1,6 +1,6 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 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, & 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 ! 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) :: PHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,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_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(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) :: 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 :: eGW(:)
double precision,allocatable :: eOld(:) double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
double precision,allocatable :: Om(:) double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:) double precision,allocatable :: XpY(:,:)
@ -100,13 +98,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Memory allocation ! 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)) 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 ! Initialization
nSCF = 0 nSCF = 0
@ -152,7 +146,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Solve the quasi-particle equation ! Solve the quasi-particle equation
eGW(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) eGW(:) = eHF(:) + SigC(:)
! Linearized or graphical solution? ! 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(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*) 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 end if

17
src/GW/evUGW.f90
View File

@ -1,6 +1,6 @@
subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & 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, & 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 ! 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) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,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) :: S(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(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 :: eOld(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: SigX(:,:)
double precision,allocatable :: SigC(:,:) double precision,allocatable :: SigC(:,:)
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) 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_bb = nS(2)
nS_sc = nS_aa + nS_bb 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), & 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)) 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 ! Initialization
nSCF = 0 nSCF = 0
@ -125,7 +116,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
Conv = 1d0 Conv = 1d0
e_diis(:,:,:) = 0d0 e_diis(:,:,:) = 0d0
error_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0
eGW(:,:) = eG0W0(:,:) eGW(:,:) = eHF(:,:)
eOld(:,:) = eGW(:,:) eOld(:,:) = eGW(:,:)
Z(:,:) = 1d0 Z(:,:) = 1d0
rcond(:) = 0d0 rcond(:) = 0d0
@ -167,7 +158,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Solve the quasi-particle equation ! ! Solve the quasi-particle equation !
!-----------------------------------! !-----------------------------------!
eGW(:,:) = eHF(:,:) + SigX(:,:) + SigC(:,:) - Vxc(:,:) eGW(:,:) = eHF(:,:) + SigC(:,:)
! Convergence criteria ! Convergence criteria

6
src/GW/unrestricted_QP_graph.f90
View File

@ -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 ! 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 integer,intent(in) :: nS
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas) 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) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin) 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) 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) 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 df = 1d0 - dsigC
w = w - f/df w = w - f/df

47
src/QuAcK/QuAcK.f90
View File

@ -28,7 +28,6 @@ program QuAcK
double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:) double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:)
double precision,allocatable :: Vxc(:,:)
logical :: doACFDT logical :: doACFDT
logical :: exchange_kernel 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), & 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), & 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 ! Read integrals
@ -217,7 +216,7 @@ program QuAcK
call wall_time(start_HF) call wall_time(start_HF)
call HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,n_diis_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, & 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) call wall_time(end_HF)
t_HF = end_HF - start_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, & 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, & 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 else
call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 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 end if
call cpu_time(end_GW) 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, & 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, & 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, & 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 else
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & 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, & 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 end if
call cpu_time(end_GW) call cpu_time(end_GW)
@ -952,15 +951,15 @@ program QuAcK
if(unrestricted) then if(unrestricted) then
call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA, & call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA, &
spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & 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, & 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) dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
else else
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 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 end if
@ -982,18 +981,18 @@ program QuAcK
if(unrestricted) then if(unrestricted) then
call evUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & call evUGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO, & 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, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, &
dipole_int_bb,PHF,cHF,epsHF,Vxc) dipole_int_bb,PHF,cHF,epsHF)
else else
call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta_GT,regGT, & 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, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO, &
PHF,cHF,epsHF,Vxc) PHF,cHF,epsHF)
end if end if
@ -1015,10 +1014,10 @@ program QuAcK
if(unrestricted) then if(unrestricted) then
call qsUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T, & 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,& 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,& 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) ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
else else
call qsGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & call qsGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
@ -1051,7 +1050,7 @@ program QuAcK
else else
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 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 end if
@ -1076,7 +1075,7 @@ program QuAcK
call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & 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, & 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 end if
call cpu_time(end_GT) call cpu_time(end_GT)