fix GW+C spectral function

src/GW/GWC_spectral_function.f90

@@ -56,8 +56,8 @@ subroutine GWC_spectral_function(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   do g=1,nGrid
     do p=nC+1,nBas-nR
 
-      ReSigC(p,g) = GW_ReSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
+      ReSigC(p,g) = GW_ReSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
-      ImSigC(p,g) = GW_ImSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
+      ImSigC(p,g) = GW_ImSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
 
     end do
   end do
@@ -91,8 +91,8 @@ subroutine GWC_spectral_function(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   do g=1,nGrid
     do p=nC+1,nBas-nR
 
-      RedSigC(p,g) = GW_RedSigC(p,0d0,eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
+      RedSigC(p,g) = GW_RedSigC(p,eHF(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
-      ImdSigC(p,g) = GW_ImdSigC(p,0d0,eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
+      ImdSigC(p,g) = GW_ImdSigC(p,eHF(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
 
     end do
   end do
@@ -106,7 +106,7 @@ subroutine GWC_spectral_function(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   do g=1,nGrid
     do p=nC+1,nBas-nR
 
-      RedSigC(p,g) = GW_RedSigC(p,eHF(p)-w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
+      RedSigC(p,g) = GW_RedSigC(p,eHF(p)+w(g),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
 
     end do
   end do

src/GW/RG0W0.f90

@@ -154,8 +154,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Cumulant expansion !
 !--------------------!
 
-  call RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
+  call RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z)
-  call GWC_spectral_function(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
+! call GWC_spectral_function(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
 
 ! Compute the RPA correlation energy
 

src/GW/RGWC.f90

@@ -1,4 +1,4 @@
-subroutine RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
+subroutine RGWC(dotest,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,Z)
 
 ! Perform GW+C calculation
 
@@ -12,11 +12,13 @@ subroutine RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
   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)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: eGW(nBas)
   double precision,intent(in)   :: Z(nBas)
 
@@ -34,6 +36,16 @@ subroutine RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
   double precision,allocatable  :: ZSat(:,:)
   double precision,allocatable  :: eSat(:,:)
 
+  double precision              :: eta
+  integer                       :: g
+  integer                       :: nGrid
+  double precision              :: wmin,wmax,dw
+  double precision,external     :: GW_ReSigC,GW_ImSigC,GW_RedSigC,GW_ImdSigC
+  double precision,allocatable  :: w(:)
+  double precision,allocatable  :: ReSigC(:,:),ImSigC(:,:)
+  double precision,allocatable  :: RedSigC(:,:),ImdSigC(:,:)
+  double precision,allocatable  :: AGWC(:,:)
+
 ! Output variables
 
 ! Hello world
@@ -165,6 +177,100 @@ subroutine RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
 
   end if
 
+! Broadening parameter
+
+  eta = 0.01d0
+
+! Construct grid
+
+  nGrid = 1000
+  allocate(w(nGrid),AGWC(nBas,nGrid))
+
+! Minimum and maximum frequency values
+
+  wmin = -5d0
+  wmax = +5d0
+  dw = (wmax - wmin)/dble(ngrid)
+
+  do g=1,nGrid
+    w(g) = wmin + dble(g)*dw
+  end do
+
+! Compute QP part of the spectral function
+
+  allocate(ReSigC(nBas,nGrid),ImSigC(nBas,nGrid))
+
+  do g=1,nGrid
+    do p=nC+1,nBas-nR
+
+      ReSigC(p,g) = GW_ReSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
+      ImSigC(p,g) = GW_ImSigC(p,eGW(p),eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
+
+    end do
+  end do
+
+  do g=1,nGrid
+    do p=nC+1,nBas-nR
+      AGWC(p,g) = ZC(p)*abs(ImSigC(p,g))/((w(g) - eGW(p))**2 + ImSigC(p,g)**2)
+    end do
+  end do
+
+  AGWC(:,:) = AGWC(:,:)/pi
+
+! Dump quantities in files as a function of w
+
+  open(unit=11 ,file='GWC_AQP.dat')
+
+  do g=1,nGrid
+    write(11,*) w(g)*HaToeV,(AGWC(p,g),p=nC+1,nBas-nR)
+  end do
+
+! Closing files
+
+  close(unit=11)
+
+! Compute cumulant part of the spectral function 
+
+  AGWC(:,:) = 0d0
+
+  if(do_hole_branch) then
+
+    do g=1,nGrid
+      do i=nC+1,nO
+        do m=1,nS
+          AGWC(i,g) = AGWC(i,g) + ZSat(i,m)*abs(ImSigC(i,g))/((w(g) - eSat(i,m))**2 + ImSigC(i,g)**2)
+        end do
+      end do
+    end do
+
+  end if
+
+  if(do_electron_branch) then
+
+    do g=1,nGrid
+      do a=nO+1,nBas-nR
+        do m=1,nS
+          AGWC(a,g) = AGWC(a,g) + ZSat(a,m)*abs(ImSigC(a,g))/((w(g) - eSat(a,m))**2 + ImSigC(a,g)**2)
+        end do
+      end do
+    end do
+
+  end if
+
+  AGWC(:,:) = AGWC(:,:)/pi
+
+! Dump quantities in files as a function of w
+
+  open(unit=12 ,file='GWC_AC.dat')
+
+  do g=1,nGrid
+    write(12,*) w(g)*HaToeV,(AGWC(p,g),p=nC+1,nBas-nR)
+  end do
+
+! Closing files
+
+  close(unit=12)
+
 ! Testing zone
 
 ! if(dotest) then