10
1
mirror of https://github.com/pfloos/quack synced 2024-11-03 12:43:48 +01:00

GW spectal function

This commit is contained in:
Pierre-Francois Loos 2024-01-29 15:23:34 +01:00
parent 7f4268a185
commit c2fc694fdb

View File

@ -24,9 +24,9 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
integer :: p,g integer :: p,g
integer :: nGrid integer :: nGrid
double precision :: wmin,wmax,dw double precision :: wmin,wmax,dw
double precision,external :: GW_SigC,GW_dSigC double precision,external :: GW_SigC,GW_ImSigC,GW_dSigC
double precision,allocatable :: w(:) double precision,allocatable :: w(:)
double precision,allocatable :: SigC(:,:) double precision,allocatable :: ReSigC(:,:),ImSigC(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
double precision,allocatable :: S(:,:) double precision,allocatable :: S(:,:)
@ -37,12 +37,13 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
! Construct grid ! Construct grid
nGrid = 5000 nGrid = 5000
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid)) allocate(w(nGrid),ReSigC(nBas,nGrid),ImSigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
! Initialize ! Initialize
SigC(:,:) = 0d0 ReSigC(:,:) = 0d0
Z(:,:) = 0d0 ImSigC(:,:) = 0d0
Z(:,:) = 0d0
! Minimum and maximum frequency values ! Minimum and maximum frequency values
@ -59,8 +60,9 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
do g=1,nGrid do g=1,nGrid
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
SigC(p,g) = GW_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) ReSigC(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) ImSigC(p,g) = GW_ImSigC(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
end do end do
@ -71,7 +73,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
do g=1,nGrid do g=1,nGrid
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2) S(p,g) = abs(ImSigC(p,g))/((w(g) - eHF(p) - ReSigC(p,g))**2 + ImSigC(p,g)**2)
end do end do
end do end do
@ -85,7 +87,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
open(unit=11 ,file='GW_A.dat') open(unit=11 ,file='GW_A.dat')
do g=1,nGrid do g=1,nGrid
write(8 ,*) w(g)*HaToeV,(SigC(p,g)*HaToeV,p=nC+1,nBas-nR) write(8 ,*) w(g)*HaToeV,(ReSigC(p,g)*HaToeV,p=nC+1,nBas-nR)
write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*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(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) write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR)