LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-02 03:15:31 +02:00
Code Releases Activity

clean up GTeh

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-21 20:38:49 +02:00
parent c476b258ff
commit 34fb993302
5 changed files with 27 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -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

6
input/options
View File

@ -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

100
src/GT/G0T0eh.f90
View File

@ -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

21
src/GT/GTeh_excitation_density.f90
View File

@ -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

8
src/GT/GTeh_self_energy_diag.f90
View File

@ -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