selfen_phon.f90 17.2 KB
Newer Older
sponce's avatar
sponce committed
1 2 3 4 5 6 7 8 9
  !                                                                            
  ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
  ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino  
  !                                                                            
  ! This file is distributed under the terms of the GNU General Public         
  ! License. See the file `LICENSE' in the root directory of the               
  ! present distribution, or http://www.gnu.org/copyleft.gpl.txt .             
  !                                                                            
  !-----------------------------------------------------------------------
10
  SUBROUTINE selfen_phon_q ( iqq, iq, totq )
sponce's avatar
sponce committed
11
  !-----------------------------------------------------------------------
sponce's avatar
sponce committed
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
  !!
  !!  compute the imaginary part of the phonon self energy due to electron-
  !!  phonon interaction in the Migdal approximation. This corresponds to 
  !!  the phonon linewidth (half width). The phonon frequency is taken into
  !!  account in the energy selection rule.
  !!
  !!  Use matrix elements, electronic eigenvalues and phonon frequencies
  !!  from ep-wannier interpolation.  This routine is similar to the one above
  !!  but it is ONLY called from within ephwann_shuffle and calculates 
  !!  the selfenergy for one phonon at a time.  Much smaller footprint on
  !!  the disk
  !!
  !!  RM 24/02/2014
  !!  redefined the size of coskkq, vkk, vkq within the fermi windwow
  !!  cleaned up the subroutine
  !!
sponce's avatar
sponce committed
28 29 30 31
  !-----------------------------------------------------------------------
  USE kinds,      ONLY : DP
  USE io_global,  ONLY : stdout
  use phcom,      ONLY : nmodes
32 33
  use epwcom,     ONLY : nbndsub, fsthick, efermi_read, fermi_energy,  &
                         eptemp, ngaussw, degaussw, shortrange,        &
34
                         nsmear, delta_smear, eps_acustic, specfun_ph, &
35
                         delta_approx, vme
sponce's avatar
sponce committed
36
  use pwcom,      ONLY : nelec, ef, isk
37 38
  use elph2,      ONLY : epf17, ibndmax, ibndmin, etf, wkf, xqf, wqf, nkqf, &
                         nkf, wf, nkqtotf, xqf, lambda_all, lambda_v_all,   &
Samuel Poncé's avatar
Samuel Poncé committed
39
                         dmef, vmef, gamma_all,gamma_v_all, efnew
40
  USE constants_epw, ONLY : ryd2mev, ryd2ev, two, zero, pi, eps4, eps6, eps8
41
  use mp,         ONLY : mp_barrier, mp_sum
Samuel Poncé's avatar
Samuel Poncé committed
42
  use mp_global,  ONLY : inter_pool_comm
sponce's avatar
sponce committed
43 44 45
  !
  implicit none
  !
46 47
  INTEGER, INTENT (in) :: iqq
  !! Current q-point index from the selecq
48 49
  INTEGER, INTENT (in) :: iq
  !! Current q-point index
50 51
  INTEGER, INTENT (in) :: totq
  !! Total number of q-points in selecq.fmt
52 53
  ! 
  ! Local variables 
sponce's avatar
sponce committed
54
  !
55 56 57 58 59 60 61 62 63 64 65 66
  INTEGER :: ik
  !! Counter on the k-point index 
  INTEGER :: ikk
  !! k-point index
  INTEGER :: ikq
  !! q-point index 
  INTEGER :: ibnd
  !! Counter on bands
  INTEGER :: jbnd
  !! Counter on bands
  INTEGER :: imode
  !! Counter on mode
67 68
  INTEGER :: jmode
  !! Counter on mode
69 70 71 72 73
  INTEGER :: fermicount
  !! Number of states on the Fermi surface
  INTEGER :: ismear
  !! Upper bounds index after k or q paral
  !! Smearing for the Gaussian function 
74 75
  INTEGER :: n
  !! Counter on number of mode degeneracies
76 77 78 79 80 81 82 83 84
  ! 
  REAL(kind=DP) :: g2
  !! Electron-phonon matrix elements squared in Ry^2
  REAL(kind=DP) :: ekk
  !! Eigen energy on the fine grid relative to the Fermi level
  REAL(kind=DP) :: ekq
  !! Eigen energy of k+q on the fine grid relative to the Fermi level
  REAL(kind=DP) :: wq
  !! Phonon frequency on the fine grid
85 86
  REAL(kind=DP) :: wq_tmp
  !! Temporary Phonon frequency on the fine grid
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
  REAL(kind=DP) :: ef0
  !! Fermi energy level
  REAL(kind=DP) :: wgkq
  !! Fermi-Dirac occupation factor $f_{nk+q}(T)$
  REAL(kind=DP) :: weight
  !! Imaginary part of the phonhon self-energy factor 
  !!$$ \pi N_q \Im(\frac{f_{nk}(T) - f_{mk+q(T)}}{\varepsilon_{nk}-\varepsilon_{mk+q}-\omega_{q\nu}+i\delta}) $$
  !! In practice the imaginary is performed with a delta Dirac
  REAL(kind=DP) :: dosef
  !! Density of state N(Ef)
  REAL(kind=DP) :: w0g1
  !! Dirac delta for the imaginary part of $\Sigma$
  REAL(kind=DP) :: w0g2
  !! Dirac delta for the imaginary part of $\Sigma$
  REAL(kind=DP) :: inv_wq
  !! $frac{1}{2\omega_{q\nu}}$ defined for efficiency reasons
  REAL(kind=DP) :: inv_eptemp0
  !! Inverse of temperature define for efficiency reasons
  REAL(kind=DP) :: g2_tmp
  !! If the phonon frequency is too small discart g
  REAL(kind=DP) :: gamma(nmodes)
  !! Gamma is the imaginary part of the phonon self-energy 
  REAL(kind=DP) :: gamma_v(nmodes)
  !! Gamma is the imaginary part of the phonon self-energy multiplied by (1-coskkq)
  REAL(kind=DP) :: coskkq(ibndmax-ibndmin+1, ibndmax-ibndmin+1)
  !! $$(v_k \cdot v_{k+q}) / |v_k|^2$$
  REAL(kind=DP) :: DDOT
  !! Dot product function
  REAL(kind=DP) :: degaussw0
  !! degaussw0 = (ismear-1) * delta_smear + degaussw
  REAL(kind=DP) :: inv_degaussw0
  !! Inverse degaussw0 for efficiency reasons
  REAL(kind=DP) :: lambda_tot
  !! Integrated lambda function
  REAL(kind=DP) :: lambda_tr_tot
  !! Integrated transport lambda function
  REAL(kind=DP) :: wgkk
  !! Fermi-Dirac occupation factor $f_{nk}(T)$
  REAL(kind=DP) :: eptemp0
  !!eptemp0   = (ismear-1) * delta_smear + eptem
  REAL(kind=DP) :: vkk(3,ibndmax-ibndmin+1)
  !! Electronic velocity $v_{nk}$
  REAL(kind=DP) :: vkq(3,ibndmax-ibndmin+1)
  !! Electronic velocity $v_{nk+q}$
131 132 133 134
  REAL(kind=DP) :: tmp
  !! Temporary value of lambda for av.
  REAL(kind=DP) :: tmp2
  !! Temporary value of lambda_v for av.
135 136 137 138
  REAL(kind=DP) :: tmp3
  !! Temporary value of lambda_v for av.
  REAL(kind=DP) :: tmp4
  !! Temporary value of lambda_v for av.
139 140 141 142
  REAL(kind=DP) :: lambda_tmp(nmodes)
  !! Temporary value of lambda for av.  
  REAL(kind=DP) :: lambda_v_tmp(nmodes)
  !! Temporary value of lambda v for av.  
143 144 145 146
  REAL(kind=DP) :: gamma_tmp(nmodes)
  !! Temporary value of gamma for av.  
  REAL(kind=DP) :: gamma_v_tmp(nmodes)
  !! Temporary value of gamma v for av.  
147 148 149 150 151 152 153 154 155 156
  REAL(kind=DP), external :: dos_ef
  !! Function to compute the Density of States at the Fermi level
  REAL(kind=DP), external :: wgauss
  !! Fermi-Dirac distribution function (when -99)
  REAL(kind=DP), external :: w0gauss
  !! This function computes the derivative of the Fermi-Dirac function
  !! It is therefore an approximation for a delta function
  REAL(kind=DP), external :: efermig
  !! Return the fermi energy
  !  
157
  IF ( iq == 1 ) THEN 
158 159 160 161 162 163 164 165 166 167
    WRITE(stdout,'(/5x,a)') repeat('=',67)
    WRITE(stdout,'(5x,"Phonon (Imaginary) Self-Energy in the Migdal Approximation")') 
    WRITE(stdout,'(5x,a/)') repeat('=',67)
    !
    IF ( fsthick.lt.1.d3 ) &
         WRITE(stdout, '(/5x,a,f10.6,a)' ) &
         'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
    WRITE(stdout, '(/5x,a,f10.6,a)' ) &
         'Golden Rule strictly enforced with T = ',eptemp * ryd2ev, ' eV'
    !
168 169
    IF ( .not. ALLOCATED (lambda_all) )   ALLOCATE( lambda_all  (nmodes, totq, nsmear) )
    IF ( .not. ALLOCATED (lambda_v_all) ) ALLOCATE( lambda_v_all(nmodes, totq, nsmear) )
170 171
    lambda_all(:,:,:)   = zero
    lambda_v_all(:,:,:) = zero
172 173
    IF ( .not. ALLOCATED (gamma_all) )    ALLOCATE( gamma_all  (nmodes, totq, nsmear) )
    IF ( .not. ALLOCATED (gamma_v_all) )  ALLOCATE( gamma_v_all(nmodes, totq, nsmear) )
174 175 176
    gamma_all(:,:,:)   = zero
    gamma_v_all(:,:,:) = zero
    !
sponce's avatar
sponce committed
177 178 179
  ENDIF
  !
  DO ismear = 1, nsmear
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
    !
    degaussw0 = (ismear-1) * delta_smear + degaussw
    eptemp0   = (ismear-1) * delta_smear + eptemp
    ! 
    ! SP: Multiplication is faster than division ==> Important if called a lot
    !     in inner loops
    inv_degaussw0 = 1.0/degaussw0
    inv_eptemp0   = 1.0/eptemp0
    !
    ! Fermi level and corresponding DOS
    !
    IF ( efermi_read ) THEN
      !
      ef0 = fermi_energy
      !
    ELSE IF (nsmear > 1) THEN
      !
      ef0 = efermig(etf,nbndsub,nkqf,nelec,wkf,degaussw0,ngaussw,0,isk)
      ! if some bands are skipped (nbndskip.neq.0), nelec has already been
      ! recalculated 
      ! in ephwann_shuffle
      !
    ELSE !SP: This is added for efficiency reason because the efermig routine is slow
      ef0 = efnew
    ENDIF
    !
    dosef = dos_ef (ngaussw, degaussw0, ef0, etf, wkf, nkqf, nbndsub)
    !  N(Ef) in the equation for lambda is the DOS per spin
    dosef = dosef / two
    !
210
    IF ( iq == 1 ) THEN 
211 212 213 214 215 216 217 218
      WRITE (stdout, 100) degaussw0 * ryd2ev, ngaussw
      WRITE (stdout, 101) dosef / ryd2ev, ef0 * ryd2ev
      !WRITE (stdout, 101) dosef / ryd2ev, ef  * ryd2ev
    ENDIF
    !
    CALL start_clock('PH SELF-ENERGY')
    !
    fermicount = 0
Samuel Poncé's avatar
Samuel Poncé committed
219 220
    wgkk = 0.0_DP
    w0g1 = 0.0_DP
221 222 223 224 225 226 227 228
    gamma(:)   = zero
    gamma_v(:) = zero
    !
    DO ik = 1, nkf
      !
      ikk = 2 * ik - 1
      ikq = ikk + 1
      ! 
Samuel Poncé's avatar
Samuel Poncé committed
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
      coskkq = zero
      ! coskkq = (vk dot vkq) / |vk|^2  appears in Grimvall 8.20
      ! this is different from :   coskkq = (vk dot vkq) / |vk||vkq|
      ! In principle the only coskkq contributing to lambda_tr are both near the
      ! Fermi surface and the magnitudes will not differ greatly between vk and vkq
      ! we may implement the approximation to the angle between k and k+q
      ! vectors also listed in Grimvall
      !
      IF (vme ) THEN 
        DO ibnd = 1, ibndmax-ibndmin+1
          DO jbnd = 1, ibndmax-ibndmin+1
            !
            ! vmef is in units of Ryd * bohr
            !
            vkk(:, ibnd ) = REAL (vmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk ) )
            vkq(:, jbnd ) = REAL (vmef (:, ibndmin-1+jbnd, ibndmin-1+jbnd, ikq ) )
            IF ( abs ( vkk(1,ibnd)**2 + vkk(2,ibnd)**2 + vkk(3,ibnd)**2) > eps4) &
                coskkq(ibnd, jbnd ) = DDOT(3, vkk(:,ibnd ), 1, vkq(:,jbnd),1)  / &
                DDOT(3, vkk(:,ibnd), 1, vkk(:,ibnd),1)
          ENDDO
sponce's avatar
sponce committed
249
        ENDDO
Samuel Poncé's avatar
Samuel Poncé committed
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
      ELSE
        DO ibnd = 1, ibndmax-ibndmin+1
          DO jbnd = 1, ibndmax-ibndmin+1
            !
            ! v_(k,i) = 1/m <ki|p|ki> = 2 * dmef (:, i,i,k)
            ! 1/m  = 2 in Rydberg atomic units
            !
            vkk(:, ibnd ) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk ) )
            vkq(:, jbnd ) = 2.0 * REAL (dmef (:, ibndmin-1+jbnd, ibndmin-1+jbnd, ikq ) )
            IF ( abs ( vkk(1,ibnd)**2 + vkk(2,ibnd)**2 + vkk(3,ibnd)**2) > eps4) &
                coskkq(ibnd, jbnd ) = DDOT(3, vkk(:,ibnd ), 1, vkq(:,jbnd),1)  / &
                DDOT(3, vkk(:,ibnd), 1, vkk(:,ibnd),1)
          ENDDO
        ENDDO
      ENDIF
265 266 267 268 269 270 271 272 273
      !
      !DBSP
      !if (ik==3) THEN
      !  print*,'vkk(:, 2)',vkk(:, 2)
      !  print*,'vkq(:, 2)',vkq(:, 2)
      !ENDIF         
      !
      ! here we must have ef, not ef0, to be consistent with ephwann_shuffle
      IF ( ( minval ( abs(etf (:, ikk) - ef) ) .lt. fsthick ) .AND. &
Samuel Poncé's avatar
Samuel Poncé committed
274
           ( minval ( abs(etf (:, ikq) - ef) ) .lt. fsthick ) ) THEN
sponce's avatar
sponce committed
275
        !
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
        fermicount = fermicount + 1
        DO imode = 1, nmodes
          !
          ! the phonon frequency
          wq = wf (imode, iq)
          !
          ! SP : We should avoid branching statements (if statements) in
          !      innerloops. Therefore we do it here.
          inv_wq =  1.0/(two * wq)
          ! the coupling from Gamma acoustic phonons is negligible
          IF ( wq .gt. eps_acustic ) THEN
            g2_tmp = 1.0
          ELSE
            g2_tmp = 0.0
          ENDIF   
          !
          DO ibnd = 1, ibndmax-ibndmin+1
            !
            !  the fermi occupation for k
            ekk = etf (ibndmin-1+ibnd, ikk) - ef0
            IF (delta_approx) THEN
              w0g1 = w0gauss ( ekk / degaussw0, 0) / degaussw0
            ELSE
              wgkk = wgauss( -ekk*inv_eptemp0, -99)
            ENDIF
            !
            DO jbnd = 1, ibndmax-ibndmin+1
sponce's avatar
sponce committed
303
              !
304 305
              !  the fermi occupation for k+q
              ekq = etf (ibndmin-1+jbnd, ikq) - ef0
sponce's avatar
sponce committed
306
              !
307 308 309
              ! here we take into account the zero-point sqrt(hbar/2M\omega)
              ! with hbar = 1 and M already contained in the eigenmodes
              ! g2 is Ry^2, wkf must already account for the spin factor
sponce's avatar
sponce committed
310
              !
311 312
              IF ( shortrange .AND. ( abs(xqf (1, iq))> eps8 .OR. abs(xqf (2, iq))> eps8 &
                 .OR. abs(xqf (3, iq))> eps8 )) THEN              
313
                ! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary 
314
                !     number, in which case its square will be a negative number. 
Samuel Poncé's avatar
Samuel Poncé committed
315
                g2 = REAL( (epf17 (jbnd, ibnd, imode, ik)**two)*inv_wq*g2_tmp ) 
sponce's avatar
sponce committed
316
              ELSE
317
                g2 = (abs(epf17 (jbnd, ibnd, imode, ik))**two)*inv_wq*g2_tmp
sponce's avatar
sponce committed
318 319
              ENDIF
              !
320 321 322 323 324 325 326 327 328 329 330 331
              IF (delta_approx) THEN 
                !
                w0g2 = w0gauss ( ekq / degaussw0, 0) / degaussw0
                ! the expression below is positive-definite, but also an
                ! approximation which neglects some fine features
                weight = pi * wq * wkf (ikk) * w0g1 * w0g2
                !
              ELSE
                !
                wgkq = wgauss( -ekq*inv_eptemp0, -99)
                !
                ! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q + id]
332
                ! This is the imaginary part of minus the phonon self-energy, sans
333 334 335 336 337
                ! the matrix elements
                !
                !weight = wkf (ikk) * (wgkk - wgkq) * &
                !     aimag ( cone / ( ekq - ekk - wq - ci * degaussw0 ) )
                !
338
                ! SP: The expression below (minus phonon self-energy) corresponds to
339 340 341 342 343
                !  = pi*k-point weight*[f(E_k) - f(E_k+q)]*delta[E_k+q - E_k - w_q]
                weight = pi * wkf (ikk) * (wgkk - wgkq)* &
                     w0gauss ( (ekq - ekk - wq) / degaussw0, 0) / degaussw0
                !
              ENDIF  
sponce's avatar
sponce committed
344
              !
345 346 347 348 349 350 351 352
              gamma   (imode) = gamma   (imode) + weight * g2 
              gamma_v (imode) = gamma_v (imode) + weight * g2 * (1-coskkq(ibnd, jbnd) ) 
              !
            ENDDO ! jbnd
            !
          ENDDO   ! ibnd
          !
        ENDDO ! loop on q-modes
sponce's avatar
sponce committed
353
        !
354 355 356 357 358 359 360 361 362 363 364 365 366
      ENDIF ! endif fsthick
      !
    ENDDO ! loop on k
    !
    CALL stop_clock('PH SELF-ENERGY')
    !
    ! collect contributions from all pools (sum over k-points)
    ! this finishes the integral over the BZ  (k)
    !
    CALL mp_sum(gamma,inter_pool_comm) 
    CALL mp_sum(gamma_v,inter_pool_comm) 
    CALL mp_sum(fermicount, inter_pool_comm)
    CALL mp_barrier(inter_pool_comm)
367 368 369 370 371 372
    ! 
    ! An average over degenerate phonon-mode is performed. 
    DO imode = 1, nmodes
      n = 0
      tmp = 0.0_DP
      tmp2 = 0.0_DP
373 374
      tmp3 = 0.0_DP
      tmp4 = 0.0_DP
375 376 377 378 379
      wq = wf (imode, iq)
      DO jmode = 1, nmodes
        wq_tmp = wf (jmode, iq)
        IF ( ABS(wq - wq_tmp) < eps6 ) THEN
          n = n + 1
380 381 382 383 384 385
          IF ( wq_tmp .gt. eps_acustic ) THEN 
            tmp  =  tmp  + gamma  ( jmode ) / pi / wq**two / dosef
            tmp2 =  tmp2 + gamma_v( jmode ) / pi / wq**two / dosef
          ENDIF
          tmp3 =  tmp3 + gamma(jmode)
          tmp4 =  tmp4 + gamma_v(jmode)
386 387 388 389
        ENDIF
      ENDDO ! jbnd
      lambda_tmp(imode)   = tmp / float(n)
      lambda_v_tmp(imode) = tmp2 / float(n)
390 391
      gamma_tmp(imode)    = tmp3 / float(n)
      gamma_v_tmp(imode)  = tmp4 / float(n)
392
    ENDDO
393
    lambda_all( :, iq, ismear )   = lambda_tmp(:)
394
    lambda_v_all( :, iq, ismear ) = lambda_v_tmp(:)
395 396 397 398
    gamma_all( :, iq, ismear )    = gamma_tmp(:)
    gamma_v_all( :, iq, ismear )  = gamma_v_tmp(:)
    lambda_tot = sum(lambda_all(:,iq,ismear))
    lambda_tr_tot = sum(lambda_v_all(:,iq,ismear))
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
    !
    WRITE(stdout,'(/5x,"ismear = ",i5," iq = ",i7," coord.: ", 3f9.5, " wt: ", f9.5)') ismear, iq, xqf(:,iq), wqf(iq)
    WRITE(stdout,'(5x,a)') repeat('-',67)
    !
    DO imode = 1, nmodes
      ! 
      wq = wf (imode, iq)
      WRITE(stdout, 102) imode, lambda_all(imode,iq,ismear),ryd2mev*gamma_all(imode,iq,ismear), ryd2mev*wq
      WRITE(stdout, 104) imode, lambda_v_all(imode,iq,ismear),ryd2mev*gamma_v_all(imode,iq,ismear), ryd2mev*wq
      !
    ENDDO
    !
    WRITE(stdout, 103) lambda_tot
    WRITE(stdout, 105) lambda_tr_tot
    ! 
414
    IF (.NOT. specfun_ph) THEN
415 416 417 418
      WRITE(stdout,'(5x,a/)') repeat('-',67)
      WRITE( stdout, '(/5x,a,i8,a,i8/)' ) &
          'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nkqtotf/2
    ENDIF
419
    !
sponce's avatar
sponce committed
420 421 422 423
  ENDDO !smears
  !
100 FORMAT(5x,'Gaussian Broadening: ',f10.6,' eV, ngauss=',i4)
101 FORMAT(5x,'DOS =',f10.6,' states/spin/eV/Unit Cell at Ef=',f10.6,' eV')
424 425
102 FORMAT(5x,'lambda___( ',i3,' )=',f15.6,'   gamma___=',f15.6,' meV','   omega=',f12.4,' meV')
103 FORMAT(5x,'lambda___( tot )=',f15.6)
sponce's avatar
sponce committed
426 427
104 FORMAT(5x,'lambda_tr( ',i3,' )=',f15.6,'   gamma_tr=',f15.6,' meV','   omega=',f12.4,' meV')
105 FORMAT(5x,'lambda_tr( tot )=',f15.6)
sponce's avatar
sponce committed
428 429 430 431
  !
  RETURN
  !
END SUBROUTINE selfen_phon_q
432
! 
sponce's avatar
sponce committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
!-----------------------------------------------------------------------
FUNCTION dos_ef_seq (ngauss, degauss, ef, et, wk, nks, nbnd)
  !-----------------------------------------------------------------------
  !
  USE kinds, ONLY : DP
  USE mp,        ONLY : mp_sum
  IMPLICIT NONE
  REAL(DP) :: dos_ef_seq
  INTEGER :: ngauss, nbnd, nks
  REAL(DP) :: et (nbnd, nks), wk (nks), ef, degauss
  !
  INTEGER :: ik, ibnd
  REAL(DP), EXTERNAL :: w0gauss
  !
  !     Compute DOS at E_F (states per Ry per unit cell)
  !
  dos_ef_seq = 0.0d0
  DO ik = 1, nks
     DO ibnd = 1, nbnd
        dos_ef_seq = dos_ef_seq + wk (ik) * w0gauss ( (et (ibnd, ik) - ef) &
             / degauss, ngauss) / degauss
     ENDDO
  ENDDO
  !
  RETURN
END FUNCTION dos_ef_seq