PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

plotting routines

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-08-24 00:02:23 +02:00
parent b96f11ad2d
commit 7efca56054
9 changed files with 221 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3 # G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F T F F F F F
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
F F F F F T F F F T F F
# * unrestricted version available # * unrestricted version available

2
input/options
View File

@ -11,7 +11,7 @@
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg
256 0.00001 T 5 F 0.0 F F 256 0.00001 T 5 F 0.0 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 F 0.1 F F 256 0.00001 T 5 F 0.0 F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA

2
src/GT/G0T0eh.f90
View File

@ -149,6 +149,8 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
end if end if
call GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR)
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies ! 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,Aph) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)

11
src/GT/GTeh_QP_graph.f90
View File

@ -1,9 +1,12 @@
subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z) subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
! Compute the graphical solution of the QP equation
implicit none implicit none
include 'parameters.h' include 'parameters.h'
! Iput variables ! Input variables
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
@ -20,6 +23,7 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
double precision,intent(in) :: eGTlin(nBas) double precision,intent(in) :: eGTlin(nBas)
! Local variables ! Local variables
integer :: p integer :: p
integer :: nIt integer :: nIt
integer,parameter :: maxIt = 64 integer,parameter :: maxIt = 64
@ -34,10 +38,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
double precision,intent(out) :: eGT(nBas) double precision,intent(out) :: eGT(nBas)
double precision,intent(out) :: Z(nBas) double precision,intent(out) :: Z(nBas)
sigC = 0d0
dsigC = 0d0
! Run Newton's algorithm to find the root ! Run Newton's algorithm to find the root
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
write(*,*) '-----------------' write(*,*) '-----------------'
@ -73,6 +75,7 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
eGT(p) = w eGT(p) = w
Z(p) = df Z(p) = df
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGT(p)*HaToeV
write(*,*) write(*,*)

98
src/GT/GTeh_plot_self_energy.f90 Normal file
View File

@ -0,0 +1,98 @@
subroutine GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR)
! Dump several GTeh quantities for external plotting
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS)
double precision,intent(in) :: rhoR(nBas,nBas,nS)
! Local variables
integer :: p,g
integer :: nGrid
double precision :: wmin,wmax,dw
double precision,external :: GTeh_SigC,GTeh_dSigC
double precision,allocatable :: w(:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: S(:,:)
! Construct grid
nGrid = 1000
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
! Initialize
SigC(:,:) = 0d0
Z(:,:) = 0d0
! Minimum and maximum frequency values
wmin = -5d0
wmax = +5d0
dw = (wmax - wmin)/dble(ngrid)
do g=1,nGrid
w(g) = wmin + dble(g)*dw
enddo
! Occupied part of the self-energy and renormalization factor
do g=1,nGrid
do p=nC+1,nBas-nR
SigC(p,g) = GTeh_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
Z(p,g) = GTeh_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
end do
end do
Z(:,:) = 1d0/(1d0 + Z(:,:))
! Compute spectral function
do g=1,nGrid
do p=nC+1,nBas-nR
S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2)
enddo
enddo
S(:,:) = S(:,:)/pi
! Dump quantities in files as a function of w
open(unit=8 ,file='GTeh_SigC.dat')
open(unit=9 ,file='GTeh_freq.dat')
open(unit=10 ,file='GTeh_Z.dat')
open(unit=11 ,file='GTeh_A.dat')
do g=1,nGrid
write(8 ,*) w(g)*HaToeV,(SigC(p,g)*HaToeV,p=nC+1,nBas-nR)
write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR)
write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR)
write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR)
enddo
! Closing files
close(unit=8)
close(unit=9)
close(unit=10)
close(unit=11)
end subroutine

6
src/GT/print_G0T0pp.f90
View File

@ -48,9 +48,9 @@ subroutine print_G0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
enddo enddo
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO energy (eV) =',eGT(HOMO)*HaToeV,' eV' write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO energy =',eGT(HOMO)*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp LUMO energy (eV) =',eGT(LUMO)*HaToeV,' eV' write(*,'(2X,A60,F15.6,A3)') 'G0T0pp LUMO energy =',eGT(LUMO)*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO-LUMO gap (eV) =',Gap*HaToeV,' eV' write(*,'(2X,A60,F15.6,A3)') 'G0T0pp HOMO-LUMO gap =',Gap*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (singlet) =',EcRPA(1),' au' write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (singlet) =',EcRPA(1),' au'
write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (triplet) =',EcRPA(2),' au' write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp correlation energy (triplet) =',EcRPA(2),' au'

2
src/GW/G0W0.f90
View File

@ -140,6 +140,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
end if end if
call GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
! Compute the RPA correlation energy ! Compute the RPA correlation energy
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)

3
src/GW/GW_QP_graph.f90
View File

@ -13,6 +13,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
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) :: Om(nS) double precision,intent(in) :: Om(nS)
@ -55,7 +56,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW,Z)
sigC = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho) sigC = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho) dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
f = w - eHF(p) - SigC f = w - eHF(p) - sigC
df = 1d0/(1d0 - dsigC) df = 1d0/(1d0 - dsigC)
w = w - df*f w = w - df*f

97
src/GW/GW_plot_self_energy.f90 Normal file
View File

@ -0,0 +1,97 @@
subroutine GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
! Dump several GW quantities for external plotting
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: p,g
integer :: nGrid
double precision :: wmin,wmax,dw
double precision,external :: GW_SigC,GW_dSigC
double precision,allocatable :: w(:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: S(:,:)
! Construct grid
nGrid = 1000
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
! Initialize
SigC(:,:) = 0d0
Z(:,:) = 0d0
! Minimum and maximum frequency values
wmin = -5d0
wmax = +5d0
dw = (wmax - wmin)/dble(ngrid)
do g=1,nGrid
w(g) = wmin + dble(g)*dw
enddo
! Occupied part of the self-energy and renormalization factor
do g=1,nGrid
do p=nC+1,nBas-nR
SigC(p,g) = GW_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
Z(p,g) = GW_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
end do
end do
Z(:,:) = 1d0/(1d0 + Z(:,:))
! Compute spectral function
do g=1,nGrid
do p=nC+1,nBas-nR
S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2)
enddo
enddo
S(:,:) = S(:,:)/pi
! Dump quantities in files as a function of w
open(unit=8 ,file='GW_SigC.dat')
open(unit=9 ,file='GW_freq.dat')
open(unit=10 ,file='GW_Z.dat')
open(unit=11 ,file='GW_A.dat')
do g=1,nGrid
write(8 ,*) w(g)*HaToeV,(SigC(p,g)*HaToeV,p=nC+1,nBas-nR)
write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR)
write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR)
write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR)
enddo
! Closing files
close(unit=8)
close(unit=9)
close(unit=10)
close(unit=11)
end subroutine