From 34fb993302a81b5b96e58e2e8542539d261780d1 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 21 Jul 2023 20:38:49 +0200 Subject: [PATCH] clean up GTeh --- input/methods | 2 +- input/options | 6 +- src/GT/G0T0eh.f90 | 100 +++-------------------------- src/GT/GTeh_excitation_density.f90 | 21 +++--- src/GT/GTeh_self_energy_diag.f90 | 8 +-- 5 files changed, 27 insertions(+), 110 deletions(-) diff --git a/input/methods b/input/methods index 47e7542..1512f13 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - T F F F F F + F F F T F T # * unrestricted version available diff --git a/input/options b/input/options index c24c4e0..232c1e1 100644 --- a/input/options +++ b/input/options @@ -5,14 +5,14 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T F T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX TDA_W reg 256 0.00001 T 5 T 0.0 F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 256 0.00001 T 5 T 0.0 F F + 256 0.00001 T 5 T 0.1 T F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA - T F F F T + F F F F T diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index b62892e..8b3a4dc 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,5 +1,5 @@ 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) ! Perform ehG0T0 calculation @@ -59,8 +59,8 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision,allocatable :: rhoL(:,:,:,:) - double precision,allocatable :: rhoR(:,:,:,:) + double precision,allocatable :: rhoL(:,:,:) + double precision,allocatable :: rhoR(:,:,:) double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) @@ -93,14 +93,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, write(*,*) end if -! Spin manifold +! Spin manifold (triplet for GTeh) ispin = 2 ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),SigX(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & - rhoL(nBas,nBas,nS,2),rhoR(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas)) + rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),eGT(nBas),eGTlin(nBas)) !-------------------! ! Compute screening ! @@ -126,9 +126,9 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, 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,Sig) -! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,Z) + + write(*,*) 'Regularization not yet implemented at the G0T0eh level!' + stop else @@ -153,14 +153,12 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, else - write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) ' *** Root search not yet implemented in G0T0eh *** ' write(*,*) - -! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGTlin,eGT) end if -! Compute the RPA correlation energy +! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) @@ -173,82 +171,4 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM) -! Deallocate memory - -! deallocate(Sig,Z,Om,XpY,XmY,rho,eGTlin) - -! Perform BSE calculation - -! if(BSE) then - -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) - -! if(exchange_kernel) then -! -! EcBSE(1) = 0.5d0*EcBSE(1) -! EcBSE(2) = 1.5d0*EcBSE(2) -! -! end if - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! Compute the BSE correlation energy via the adiabatic connection - -! if(doACFDT) then - -! write(*,*) '-------------------------------------------------------------' -! write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy ' -! write(*,*) '-------------------------------------------------------------' -! write(*,*) - -! if(doXBS) then - -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) - -! end if - -! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - -! end if - -! if(ppBSE) then - -! call Bethe_Salpeter_pp(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - -! if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eGW) -! if(BSE) call ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,Om,rho) - -! if(BSE) call XBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) - - end subroutine diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90 index 11da30e..3f3f867 100644 --- a/src/GT/GTeh_excitation_density.f90 +++ b/src/GT/GTeh_excitation_density.f90 @@ -14,17 +14,17 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) ! Local variables integer :: m,jb,p,q,j,b - double precision :: X,Y,Xt,Yt + double precision :: X,Y ! Output variables - double precision,intent(out) :: rhoL(nBas,nBas,nS,2) - double precision,intent(out) :: rhoR(nBas,nBas,nS,2) + double precision,intent(out) :: rhoL(nBas,nBas,nS) + double precision,intent(out) :: rhoR(nBas,nBas,nS) ! Initialization - rhoL(:,:,:,:) = 0d0 - rhoR(:,:,:,:) = 0d0 + rhoL(:,:,:) = 0d0 + rhoR(:,:,:) = 0d0 !$OMP PARALLEL & !$OMP SHARED(nC,nBas,nR,nO,nS,rhoL,rhoR,ERI,XpY,XmY) & @@ -42,14 +42,11 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) X = 0.5d0*(XpY(m,jb) + XmY(m,jb)) Y = 0.5d0*(XpY(m,jb) - XmY(m,jb)) - Xt = 0.5d0*(XpY(jb,m) + XmY(jb,m)) - Yt = 0.5d0*(XpY(jb,m) - XmY(jb,m)) + rhoL(p,q,m) = rhoL(p,q,m) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y +! rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,b,j,q)*X + ERI(p,j,b,q)*Y - rhoL(p,q,m,1) = rhoL(p,q,m,1) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y - rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,b,j,q)*X + ERI(p,j,b,q)*Y - - rhoR(p,q,m,1) = rhoR(p,q,m,1) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y - rhoR(p,q,m,2) = rhoR(p,q,m,2) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y + rhoR(p,q,m) = rhoR(p,q,m) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y +! rhoR(p,q,m,2) = rhoR(p,q,m,2) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y enddo enddo diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 index 41df1a1..1fcb612 100644 --- a/src/GT/GTeh_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -16,8 +16,8 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rhoL(nBas,nBas,nS,2) - double precision,intent(in) :: rhoR(nBas,nBas,nS,2) + double precision,intent(in) :: rhoL(nBas,nBas,nS) + double precision,intent(in) :: rhoR(nBas,nBas,nS) ! Local variables @@ -46,7 +46,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(i) + Om(m) - num = rhoL(p,i,m,2)*rhoR(i,p,m,1) + num = rhoL(i,p,m)*rhoR(i,p,m) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -61,7 +61,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig do m=1,nS eps = e(p) - e(a) - Om(m) - num = rhoL(p,a,m,1)*rhoR(a,p,m,2) + num = rhoL(p,a,m)*rhoR(p,a,m) Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2