PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

fix QP graph

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-02-15 22:15:30 +01:00
parent a2155a4ae6
commit ba7a2349af
4 changed files with 21 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/MBPT/G0W0.f90
View File

@ -155,7 +155,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eG0W0lin,eG0W0)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eG0W0lin,eG0W0)
! Find all the roots of the QP equation if necessary

16
src/MBPT/QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW)
! Compute the graphical solution of the QP equation
@ -15,6 +15,8 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
integer,intent(in) :: nS
double precision,intent(in) :: eta
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) :: rho(nBas,nBas,nS)
@ -27,7 +29,7 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: SigmaC,dSigmaC
double precision :: sig,dsig
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
@ -52,14 +54,14 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
nIt = nIt + 1
sig = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
dsig = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
f = w - eHF(p) - sig
df = 1d0 - dsig
sigC = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
f = w - eHF(p) - SigX(p) - sigC + Vxc(p)
df = 1d0 - dsigC
w = w - f/df
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
end do

4
src/MBPT/UG0W0.f90
View File

@ -165,8 +165,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Find graphical solution of the QP equation
do is=1,nspin
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))
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),SigX(:,is),Vxc(:,is), &
OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do
end if

16
src/MBPT/unrestricted_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW)
! Compute the graphical solution of the QP equation
@ -15,6 +15,8 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eG
integer,intent(in) :: nS
double precision,intent(in) :: eta
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) :: rho(nBas,nBas,nS,nspin)
@ -27,7 +29,7 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eG
integer,parameter :: maxIt = 10
double precision,parameter :: thresh = 1d-6
double precision,external :: USigmaC,dUSigmaC
double precision :: sig,dsig
double precision :: sigC,dsigC
double precision :: f,df
double precision :: w
@ -52,14 +54,14 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eG
nIt = nIt + 1
sig = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
dsig = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
f = w - eHF(p) - sig
df = 1d0 - dsig
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)
f = w - eHF(p) - SigX(p) + Vxc(p) - sigC
df = 1d0 - dsigC
w = w - f/df
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
end do