mirror of https://github.com/pfloos/quack
183 lines
4.7 KiB
Fortran
183 lines
4.7 KiB
Fortran
subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW)
|
|
|
|
! Compute all the roots of the QP equation for each orbital
|
|
|
|
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) :: Omega(nS)
|
|
double precision,intent(in) :: rho(nBas,nBas,nS)
|
|
|
|
! Local variables
|
|
|
|
integer :: i,j,a,b,x,jb,g
|
|
integer,parameter :: nGrid = 1000000
|
|
double precision,parameter :: wmin = -50d0
|
|
double precision,parameter :: wmax = +50d0
|
|
double precision :: dw
|
|
double precision :: eps
|
|
double precision :: left_sign,right_sign
|
|
double precision :: left_QP,right_QP
|
|
double precision,allocatable :: w(:)
|
|
double precision,allocatable :: SigC(:)
|
|
double precision,allocatable :: dSigC(:)
|
|
|
|
integer :: nIt
|
|
integer,parameter :: maxIt = 64
|
|
double precision,parameter :: thresh = 1d-6
|
|
double precision,external :: SigmaC,dSigmaC
|
|
double precision :: sig,dsig
|
|
double precision :: sat,dsat
|
|
double precision :: om
|
|
integer :: nRoot
|
|
double precision :: Z
|
|
double precision :: Ztot
|
|
double precision :: QP
|
|
|
|
! Output variables
|
|
|
|
double precision,intent(out) :: eGW(nBas)
|
|
|
|
! Construct grid
|
|
|
|
allocate(w(nGrid),SigC(nGrid),dSigC(nGrid))
|
|
|
|
! Minimum and maximum frequency values
|
|
|
|
dw = (wmax - wmin)/dble(nGrid)
|
|
|
|
do g=1,nGrid
|
|
w(g) = wmin + dble(g)*dw
|
|
enddo
|
|
|
|
! Main loop over the orbitals
|
|
|
|
do x=nC+1,nBas-nR
|
|
|
|
SigC(:) = 0d0
|
|
dSigC(:) = 0d0
|
|
|
|
! Loop over grid points
|
|
|
|
do g=1,nGrid
|
|
|
|
! Occupied part of the self-energy and spectral weight
|
|
|
|
do i=nC+1,nO
|
|
jb = 0
|
|
do j=nC+1,nO
|
|
do b=nO+1,nBas-nR
|
|
|
|
jb = jb + 1
|
|
eps = w(g) - eHF(i) + Omega(jb)
|
|
SigC(g) = SigC(g) + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
|
|
dSigC(g) = dSigC(g) - 2d0*rho(x,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
! Virtual part of the self-energy and spectral weight
|
|
|
|
do a=nO+1,nBas-nR
|
|
jb = 0
|
|
do j=nC+1,nO
|
|
do b=nO+1,nBas-nR
|
|
|
|
jb = jb + 1
|
|
eps = w(g) - eHF(a) - Omega(jb)
|
|
SigC(g) = SigC(g) + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
|
|
dSigC(g) = dSigC(g) - 2d0*rho(x,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
enddo
|
|
|
|
! Find the zeros of the QP equation
|
|
|
|
nRoot = 0
|
|
QP = 0d0
|
|
Ztot = 0d0
|
|
|
|
write(*,*) '-----------------'
|
|
write(*,'(A10,I3)') 'Orbital ',x
|
|
write(*,*) '-----------------'
|
|
write(*,*)
|
|
left_QP = w(1) - eHF(x) - SigC(1)
|
|
left_sign = sign(1d0,left_QP)
|
|
|
|
do g=2,nGrid
|
|
|
|
right_QP = w(g) - eHF(x) - SigC(g)
|
|
right_sign = sign(1d0,right_QP)
|
|
|
|
if(left_sign /= right_sign .and. left_sign == sign(1d0,-1d0)) then
|
|
|
|
nRoot = nRoot + 1
|
|
write(*,'(A20,I6,F10.6,F10.6,I6,F10.6,F10.6,F10.6)') &
|
|
'root right here!',g-1,w(g-1),left_QP,g,w(g),right_QP,1d0/(1d0-dSigC(g))
|
|
|
|
! Run Newton's algorithm to find the root
|
|
|
|
om = w(g-1)
|
|
nIt = 0
|
|
sat = 1d0
|
|
|
|
do while (abs(sat) > thresh .and. nIt < maxIt)
|
|
|
|
nIt = nIt + 1
|
|
|
|
sig = SigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
|
|
dsig = dSigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
|
|
sat = om - eHF(x) - sig
|
|
dsat = 1d0 - dsig
|
|
Z = 1d0/dsat
|
|
|
|
write(*,'(A3,I3,A1,1X,4F15.9)') 'It.',nIt,':',om,sat,sig,Z
|
|
|
|
om = om - sat/dsat
|
|
|
|
end do
|
|
|
|
if(nIt == maxIt) then
|
|
|
|
write(*,*) 'Newton root search has not converged!'
|
|
|
|
else
|
|
|
|
Ztot = Ztot + Z
|
|
QP = QP + Z*om
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
left_sign = right_sign
|
|
left_QP = right_QP
|
|
|
|
end do
|
|
|
|
eGW(x) = QP/Ztot
|
|
|
|
write(*,*)
|
|
write(*,'(A32,I3,A1,I3)') 'Number of roots for orbital ',x,':',nRoot
|
|
write(*,'(A32,F16.10)') 'Total weight ',Ztot
|
|
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(x)*HaToeV
|
|
write(*,*)
|
|
|
|
end do
|
|
|
|
end subroutine QP_roots
|