ephwann_shuffle.f90 84.9 KB
Newer Older
sponce's avatar
sponce committed
1 2 3 4 5 6 7 8 9 10 11
  !                                                                            
  ! 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 .             
  !                                                                            
  !----------------------------------------------------------------------
  SUBROUTINE ephwann_shuffle( nqc, xqc )
  !---------------------------------------------------------------------
12 13 14 15 16 17 18 19 20 21 22 23 24
  !!
  !!  Wannier interpolation of electron-phonon vertex
  !!
  !!  Scalar implementation   Feb 2006
  !!  Parallel version        May 2006
  !!  Disentenglement         Oct 2006
  !!  Compact formalism       Dec 2006
  !!  Phonon irreducible zone Mar 2007
  !!
  !!  No ultrasoft now
  !!  No spin polarization
  !!
  !!  RM - add noncolin case
sponce's avatar
sponce committed
25 26
  !-----------------------------------------------------------------------
  !
Samuel Poncé's avatar
Samuel Poncé committed
27
  USE kinds,         ONLY : DP, i4b
28
  USE pwcom,         ONLY : nbnd, nks, nkstot, isk, et, xk, ef,  nelec
29
  USE cell_base,     ONLY : at, bg, omega, alat
sponce's avatar
sponce committed
30
  USE start_k,       ONLY : nk1, nk2, nk3
31
  USE ions_base,     ONLY : nat, amass, ityp, tau
32
  USE phcom,         ONLY : nq1, nq2, nq3, nmodes
33
  USE epwcom,        ONLY : nbndsub, fsthick, epwread, longrange,               &
34
                            epwwrite, ngaussw, degaussw, lpolar, lifc, lscreen, &
35
                            nbndskip, etf_mem, scr_typ,                         &
36
                            elecselfen, phonselfen, nest_fn, a2f, specfun_ph,   &
37
                            vme, eig_read, ephwrite, nkf1, nkf2, nkf3,          & 
38
                            efermi_read, fermi_energy, specfun_el, band_plot,   &
39 40 41
                            scattering, nstemp, int_mob, scissor, carrier,      &
                            iterative_bte, longrange, scatread, nqf1, prtgkk,   &
                            nqf2, nqf3, mp_mesh_k, restart, ncarrier, plselfen, &
42 43
                            specfun_pl, lindabs, mob_maxiter, use_ws,           &
                            epmatkqread, selecqread
sponce's avatar
sponce committed
44
  USE noncollin_module, ONLY : noncolin
45 46
  USE constants_epw, ONLY : ryd2ev, ryd2mev, one, two, zero, czero,             &
                            twopi, ci, kelvin2eV, eps6, eps8 
47
  USE io_files,      ONLY : prefix, diropn, tmp_dir
sponce's avatar
sponce committed
48
  USE io_global,     ONLY : stdout, ionode
49
  USE io_epw,        ONLY : lambda_phself, linewidth_phself, iunepmatwe,        &
50
                            iunepmatwp, crystal, iunepmatwp2, iunrestart
51 52
  USE elph2,         ONLY : cu, cuq, lwin, lwinq, map_rebal, map_rebal_inv,     &
                            chw, chw_ks, cvmew, cdmew, rdw,                     &
53 54 55 56
                            epmatwp, epmatq, wf, etf, etf_k, etf_ks, xqf, xkf,  &
                            wkf, dynq, nqtotf, nkqf, epf17, nkf, nqf, et_ks,    &
                            ibndmin, ibndmax, lambda_all, dmec, dmef, vmef,     &
                            sigmai_all, sigmai_mode, gamma_all, epsi, zstar,    &
57
                            efnew, sigmar_all, zi_all, nkqtotf, eps_rpa,        &
58
                            sigmar_all, zi_allvb, inv_tau_all,                  &
59
                            inv_tau_allcb, zi_allcb, exband
60
  USE transportcom,  ONLY : transp_temp, mobilityh_save, mobilityel_save, lower_bnd, &
61 62 63
                            upper_bnd 
  USE wan2bloch,     ONLY : dmewan2bloch, hamwan2bloch, dynwan2bloch,           &
                            ephwan2blochp, ephwan2bloch, vmewan2bloch,          &
64
                            dynifc2blochf, dynifc2blochc 
65 66
  USE bloch2wan,     ONLY : hambloch2wan, dmebloch2wan, dynbloch2wan,           &
                            vmebloch2wan, ephbloch2wane, ephbloch2wanp,         &
67
                            ephbloch2wanp_mem
68
  USE wigner,        ONLY : wigner_seitz_wrap
69
  USE io_eliashberg, ONLY : write_ephmat, count_kpoints, kmesh_fine, kqmap_fine
70
  USE transport,     ONLY : transport_coeffs, scattering_rate_q, qwindow
71
  USE printing,      ONLY : print_gkk
72 73 74 75
  USE io_scattering, ONLY : electron_read, tau_read, iter_open
  USE transport_iter,ONLY : iter_restart
  USE close_epw,     ONLY : iter_close
  USE division,      ONLY : fkbounds
sponce's avatar
sponce committed
76 77
  USE mp,            ONLY : mp_barrier, mp_bcast, mp_sum
  USE io_global,     ONLY : ionode_id
78
  USE mp_global,     ONLY : inter_pool_comm, intra_pool_comm, root_pool
79
  USE mp_world,      ONLY : mpime, world_comm
sponce's avatar
sponce committed
80
#if defined(__MPI)
81 82
  USE parallel_include, ONLY : MPI_MODE_RDONLY, MPI_INFO_NULL, MPI_OFFSET_KIND, &
                               MPI_OFFSET
sponce's avatar
sponce committed
83
#endif
sponce's avatar
sponce committed
84 85 86
  !
  implicit none
  !
87 88
  INTEGER, INTENT (in) :: nqc
  !! number of qpoints in the coarse grid
sponce's avatar
sponce committed
89
  !
90 91 92 93 94 95 96 97
  REAL(kind=DP), INTENT (in) :: xqc(3,nqc)
  !! qpoint list, coarse mesh
  ! 
  ! Local  variables
  LOGICAL :: already_skipped
  !! Skipping band during the Wannierization
  LOGICAL :: exst
  !! If the file exist
98 99 100 101
  LOGICAL :: first_cycle
  !! Check wheter this is the first cycle after a restart. 
  LOGICAL :: first_time
  !! Check wheter this is the first timeafter a restart. 
102 103
  LOGICAL :: homogeneous
  !! Check if the k and q grids are homogenous and commensurate.
104 105 106 107 108 109
  !
  CHARACTER (len=256) :: filint
  !! Name of the file to write/read 
  CHARACTER (len=30)  :: myfmt
  !! Variable used for formatting output
  ! 
sponce's avatar
sponce committed
110 111
  INTEGER :: ios
  !! integer variable for I/O control
112 113
  INTEGER :: iq 
  !! Counter on coarse q-point grid
114 115
  INTEGER :: iqq
  !! Counter on coarse q-point grid  
116 117
  INTEGER :: iq_restart
  !! Counter on coarse q-point grid
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
  INTEGER :: ik
  !! Counter on coarse k-point grid
  INTEGER :: ikk
  !! Counter on k-point when you have paired k and q
  INTEGER :: ikq
  !! Paired counter so that q is adjacent to its k
  INTEGER :: ibnd
  !! Counter on band
  INTEGER :: jbnd
  !! Counter on band
  INTEGER :: imode
  !! Counter on mode
  INTEGER :: na
  !! Counter on atom
  INTEGER :: mu
  !! counter on mode
  INTEGER :: nu
  !! counter on mode
  INTEGER :: fermicount
  !! Number of states at the Fermi level
  INTEGER :: lrepmatw
  !! record length while reading file
140 141 142 143 144 145 146 147 148 149 150 151
  INTEGER :: ikx
  !! Counter on the coase k-grid
  INTEGER :: ikfx 
  !! Counter on the fine k-grid. 
  INTEGER :: xkk1, xkq1
  !! Integer of xkk when multiplied by nkf/nk
  INTEGER :: xkk2, xkq2
  !! Integer of xkk when multiplied by nkf/nk
  INTEGER :: xkk3, xkq3
  !! Integer of xkk when multiplied by nkf/nk
  INTEGER :: ir
  !! Counter for WS loop
152 153
  INTEGER :: nrws
  !! Number of real-space Wigner-Seitz
sponce's avatar
sponce committed
154 155
  INTEGER :: valueRSS(2)
  !! Return virtual and resisdent memory from system
156 157
  INTEGER :: ierr
  !! Error status
158 159 160 161 162 163
  INTEGER :: nrr_k 
  !! Number of WS points for electrons
  INTEGER :: nrr_q
  !! Number of WS points for phonons
  INTEGER :: nrr_g
  !! Number of WS points for electron-phonons
164 165
  INTEGER :: dims
  !! Dims is either nbndsub if use_ws or 1 if not
166 167
  INTEGER :: dims2
  !! Dims is either nat if use_ws or 1 if not
168 169 170 171
  INTEGER :: iw 
  !! Counter on bands when use_ws == .true.
  INTEGER :: iw2
  !! Counter on bands when use_ws == .true.
172 173 174 175 176 177 178 179
  INTEGER :: iter
  !! Current iteration number
  INTEGER :: itemp
  !! Temperature index
  INTEGER :: icbm
  !! Index of the CBM
  INTEGER :: totq
  !! Total number of q-points within the fsthick window. 
180 181 182 183 184 185 186
  INTEGER, ALLOCATABLE :: irvec_k(:,:)
  !! integer components of the ir-th Wigner-Seitz grid point in the basis
  !! of the lattice vectors for electrons
  INTEGER, ALLOCATABLE :: irvec_q(:,:)
  !! integer components of the ir-th Wigner-Seitz grid point for phonons
  INTEGER, ALLOCATABLE :: irvec_g(:,:)
  !! integer components of the ir-th Wigner-Seitz grid point for electron-phonon
187
  INTEGER, ALLOCATABLE :: ndegen_k (:,:,:)
188 189 190 191
  !! Wigner-Seitz number of degenerescence (weights) for the electrons grid
  INTEGER, ALLOCATABLE :: ndegen_q (:,:,:)
  !! Wigner-Seitz weights for the phonon grid that depend on 
  !! atomic positions $R + \tau(nb) - \tau(na)$
192
  INTEGER, ALLOCATABLE :: ndegen_g (:,:,:,:)
193 194
  !! Wigner-Seitz weights for the electron-phonon grid that depend on 
  !! atomic positions $R - \tau(na)$
195 196
  INTEGER, ALLOCATABLE :: selecq(:)
  !! Selected q-points within the fsthick window
197 198
  INTEGER, PARAMETER :: nrwsx=200
  !! Maximum number of real-space Wigner-Seitz
199 200 201 202 203 204 205 206 207 208 209
#if defined(__MPI)
  INTEGER (kind=MPI_OFFSET_KIND) :: ind_tot
  INTEGER (kind=MPI_OFFSET_KIND) :: ind_totcb
  INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw2
  INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw4
  INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw5
  INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw6
  !! Offset to tell where to start reading the file
  INTEGER (kind=MPI_OFFSET_KIND) :: lsize
  !! Offset to tell where to start reading the file
#else
Samuel Poncé's avatar
Samuel Poncé committed
210 211 212 213 214 215
  INTEGER :: ind_tot
  INTEGER :: ind_totcb
  INTEGER :: lrepmatw2
  INTEGER :: lrepmatw4
  INTEGER :: lrepmatw5
  INTEGER :: lrepmatw6
216
  !! Offset to tell where to start reading the file
Samuel Poncé's avatar
Samuel Poncé committed
217
  INTEGER :: lsize
218
#endif
219
  !  
220 221
  REAL(kind=DP) :: rdotk_scal
  !! Real (instead of array) for $r\cdot k$
222 223 224 225 226 227 228 229
  REAL(kind=DP) :: xxq(3)
  !! Current q-point 
  REAL(kind=DP) :: xxk(3)
  !! Current k-point on the fine grid
  REAL(kind=DP) :: xkk(3)
  !! Current k-point on the fine grid
  REAL(kind=DP) :: xkq(3)
  !! Current k+q point on the fine grid
230 231 232 233
  REAL(kind=DP) :: rws(0:3,nrwsx)
  !! Real-space wigner-Seitz vectors
  REAL(kind=DP) :: atws(3,3)
  !! Maximum vector: at*nq
234 235
  REAL(kind=DP) :: w_centers(3,nbndsub)
  !! Wannier centers  
236 237 238 239 240 241 242 243 244 245
  REAL(KIND=DP) :: etemp
  !! Temperature in Ry (this includes division by kb)
  REAL(KIND=DP) :: ef0(nstemp)
  !! Fermi level for the temperature itemp  
  REAL(KIND=DP) :: efcb(nstemp)
  !! Second Fermi level for the temperature itemp  
  REAL(KIND=DP) :: dummy(3)
  !! Dummy variable
  REAL(KIND=DP), EXTERNAL :: fermicarrier
  !! Function that returns the Fermi level so that n=p (if int_mob = .true.)  
246 247 248 249
  REAL(kind=DP), EXTERNAL :: efermig
  !! External function to calculate the fermi energy
  REAL(kind=DP), EXTERNAL :: efermig_seq
  !! Same but in sequential
250 251
  REAL(kind=DP), ALLOCATABLE :: etf_all(:,:)
  !! Eigen-energies on the fine grid collected from all pools in parallel case
252 253
  REAL(kind=DP), ALLOCATABLE :: w2 (:)
  !! Interpolated phonon frequency
254 255 256 257
  REAL(kind=DP), ALLOCATABLE :: irvec_r (:,:)
  !! Wigner-Size supercell vectors, store in real instead of integer
  REAL(kind=DP), ALLOCATABLE :: rdotk(:)
  !! $r\cdot k$
258 259
  REAL(kind=DP), ALLOCATABLE :: rdotk2(:)
  !! $r\cdot k$
260 261 262 263 264 265
  REAL(kind=DP), ALLOCATABLE :: wslen_k(:)
  !! real-space length for electrons, in units of alat
  REAL(kind=DP), ALLOCATABLE :: wslen_q(:)
  !! real-space length for phonons, in units of alat
  REAL(kind=DP), ALLOCATABLE :: wslen_g(:)
  !! real-space length for electron-phonons, in units of alat
266 267 268 269
  REAL(kind=DP), ALLOCATABLE :: vkk_all(:,:,:)
  !! velocity from all the k-point
  REAL(kind=DP), ALLOCATABLE :: wkf_all(:)
  !! k-point weights for all the k-points
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286
  !
  COMPLEX(kind=DP), ALLOCATABLE :: epmatwe  (:,:,:,:,:)
  !! e-p matrix  in wannier basis - electrons
  COMPLEX(kind=DP), ALLOCATABLE :: epmatwe_mem  (:,:,:,:)
  !! e-p matrix  in wannier basis - electrons (written on disk)
  COMPLEX(kind=DP), ALLOCATABLE :: epmatwef (:,:,:,:)
  !! e-p matrix  in el wannier - fine Bloch phonon grid
  COMPLEX(kind=DP), ALLOCATABLE :: epmatf( :, :, :)
  !! e-p matrix  in smooth Bloch basis, fine mesh
  COMPLEX(kind=DP), ALLOCATABLE :: cufkk ( :, :)
  !! Rotation matrix, fine mesh, points k
  COMPLEX(kind=DP), ALLOCATABLE :: cufkq ( :, :)
  !! the same, for points k+q
  COMPLEX(kind=DP), ALLOCATABLE :: uf( :, :)
  !! Rotation matrix for phonons
  COMPLEX(kind=DP), ALLOCATABLE :: bmatf ( :, :)
  !! overlap U_k+q U_k^\dagger in smooth Bloch basis, fine mesh
287
  COMPLEX(kind=DP), ALLOCATABLE :: cfac(:,:,:)
288
  !! Used to store $e^{2\pi r \cdot k}$ exponential 
289
  COMPLEX(kind=DP), ALLOCATABLE :: cfacq(:,:,:)
290
  !! Used to store $e^{2\pi r \cdot k+q}$ exponential
sponce's avatar
sponce committed
291 292
  ! 
  IF (nbndsub.ne.nbnd) &
293
       WRITE(stdout, '(/,5x,a,i4)' ) 'Band disentanglement is used:  nbndsub = ', nbndsub
sponce's avatar
sponce committed
294 295 296 297 298
  !
  ALLOCATE ( cu ( nbnd, nbndsub, nks), & 
             cuq ( nbnd, nbndsub, nks), & 
             lwin ( nbnd, nks ), &
             lwinq ( nbnd, nks ), &
299
             exband ( nbnd ) ) 
sponce's avatar
sponce committed
300 301 302
  !
  CALL start_clock ( 'ephwann' )
  !
303
  IF ( epwread ) THEN
304 305 306 307 308
    !
    ! Might have been pre-allocate depending of the restart configuration 
    IF(ALLOCATED(tau))  DEALLOCATE( tau )
    IF(ALLOCATED(ityp)) DEALLOCATE( ityp )
    IF(ALLOCATED(w2))   DEALLOCATE( w2 )
309 310 311 312 313 314 315 316 317 318
    ! 
    ! We need some crystal info
    IF (mpime.eq.ionode_id) THEN
      !
      OPEN(unit=crystal,file='crystal.fmt',status='old',iostat=ios)
      READ (crystal,*) nat
      READ (crystal,*) nmodes
      READ (crystal,*) nelec
      READ (crystal,*) at
      READ (crystal,*) bg
319 320
      READ (crystal,*) omega
      READ (crystal,*) alat
321
      ALLOCATE( tau( 3, nat ) )
322
      READ (crystal,*) tau
323 324 325
      READ (crystal,*) amass
      ALLOCATE( ityp( nat ) )
      READ (crystal,*) ityp
sponce's avatar
sponce committed
326
      READ (crystal,*) isk
327
      READ (crystal,*) noncolin
328
      READ (crystal,*) w_centers
329 330
      ! 
    ENDIF
331 332
    CALL mp_bcast (nat     , ionode_id, inter_pool_comm)
    CALL mp_bcast (nat     , root_pool, intra_pool_comm)  
333
    IF (mpime /= ionode_id) ALLOCATE( ityp( nat ) )
334 335 336 337 338 339 340 341 342 343 344 345
    CALL mp_bcast (nmodes  , ionode_id, inter_pool_comm)
    CALL mp_bcast (nmodes  , root_pool, intra_pool_comm)  
    CALL mp_bcast (nelec   , ionode_id, inter_pool_comm)
    CALL mp_bcast (nelec   , root_pool, intra_pool_comm)  
    CALL mp_bcast (at      , ionode_id, inter_pool_comm)
    CALL mp_bcast (at      , root_pool, intra_pool_comm)  
    CALL mp_bcast (bg      , ionode_id, inter_pool_comm)
    CALL mp_bcast (bg      , root_pool, intra_pool_comm)  
    CALL mp_bcast (omega   , ionode_id, inter_pool_comm)
    CALL mp_bcast (omega   , root_pool, intra_pool_comm)  
    CALL mp_bcast (alat    , ionode_id, inter_pool_comm)
    CALL mp_bcast (alat    , root_pool, intra_pool_comm)  
346
    IF (mpime /= ionode_id) ALLOCATE( tau( 3, nat ) )
347 348 349 350 351 352 353 354
    CALL mp_bcast (tau     , ionode_id, inter_pool_comm)
    CALL mp_bcast (tau     , root_pool, intra_pool_comm)  
    CALL mp_bcast (amass   , ionode_id, inter_pool_comm)
    CALL mp_bcast (amass   , root_pool, intra_pool_comm)  
    CALL mp_bcast (ityp    , ionode_id, inter_pool_comm)
    CALL mp_bcast (ityp    , root_pool, intra_pool_comm)  
    CALL mp_bcast (isk     , ionode_id, inter_pool_comm)
    CALL mp_bcast (isk     , root_pool, intra_pool_comm)  
355
    CALL mp_bcast (noncolin, ionode_id, inter_pool_comm)
356
    CALL mp_bcast (noncolin, root_pool, intra_pool_comm)  
357 358
    CALL mp_bcast (w_centers, ionode_id, inter_pool_comm)
    CALL mp_bcast (w_centers, root_pool, intra_pool_comm)
359 360 361 362 363 364 365 366
    IF (mpime.eq.ionode_id) THEN
      CLOSE(crystal)
    ENDIF
    CALL mp_barrier(inter_pool_comm)
    ! 
  ELSE
    continue
  ENDIF
sponce's avatar
sponce committed
367
  !
368 369
  ALLOCATE( w2( 3*nat) )
  !
370
  ! Determine Wigner-Seitz points
371 372 373 374 375 376 377
  ! 
  ! For this we need the Wannier centers
  ! w_centers is allocated inside loadumat
  IF (.not. epwread) THEN
    xxq = 0.d0
    CALL loadumat( nbnd, nbndsub, nks, nkstot, xxq, cu, cuq, lwin, lwinq, exband, w_centers )
  ENDIF
sponce's avatar
sponce committed
378
  !
379 380
  ! Inside we allocate irvec_k, irvec_q, irvec_g, ndegen_k, ndegen_q, ndegen_g,
  !                    wslen_k,  wslen_q,  wslen_g  
381 382
  IF (use_ws) THEN
    ! Use Wannier-centers to contstruct the WS for electonic part and el-ph part
383 384 385
    ! Use atomic position to contstruct the WS for the phonon part
    dims  = nbndsub
    dims2 = nat
386
    CALL wigner_seitz_wrap ( nk1, nk2, nk3, nq1, nq2, nq3, irvec_k, irvec_q, irvec_g, &
387 388
                             ndegen_k, ndegen_q, ndegen_g, wslen_k, wslen_q, wslen_g, &
                             w_centers, dims, tau, dims2 )
389
  ELSE
390 391 392
    ! Center the WS at Gamma for electonic part, the phonon part and el-ph part
    dims  = 1
    dims2 = 1
393 394
    dummy(:) = (/0.0,0.0,0.0/)
    CALL wigner_seitz_wrap ( nk1, nk2, nk3, nq1, nq2, nq3, irvec_k, irvec_q, irvec_g, &
395 396
                             ndegen_k, ndegen_q, ndegen_g, wslen_k, wslen_q, wslen_g, &
                             dummy, dims, dummy, dims2 )
397
  ENDIF
398
  ! 
399 400 401 402
  ! Determine the size of the respective WS sets based on the length of the matrices
  nrr_k = SIZE(irvec_k(1,:))
  nrr_q = SIZE(irvec_q(1,:))
  nrr_g = SIZE(irvec_g(1,:))
403 404 405 406 407
  IF (use_ws) THEN 
    WRITE(stdout, '(5x,a)' )    'Construct the Wigner-Seitz cell using Wannier centers and atomic positions '
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for electrons ',nrr_k
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for phonons ',nrr_q
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for electron-phonon ',nrr_g
Samuel Poncé's avatar
Samuel Poncé committed
408
    WRITE(stdout, '(5x,a,i8)' ) 'Maximum number of cores for efficient parallelization ',nrr_g * nat
409 410 411 412 413
  ELSE
    WRITE(stdout, '(5x,a)' )    'Use zone-centred Wigner-Seitz cells '
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for electrons ',nrr_k
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for phonons ',nrr_q
    WRITE(stdout, '(5x,a,i8)' ) 'Number of WS vectors for electron-phonon ',nrr_g
Samuel Poncé's avatar
Samuel Poncé committed
414
    WRITE(stdout, '(5x,a,i8)' ) 'Maximum number of cores for efficient parallelization ',nrr_g * nmodes
415 416
    WRITE(stdout, '(5x,a)' )    'Results may improve by using use_ws == .true. '
  ENDIF
sponce's avatar
sponce committed
417
  !
418 419
#ifndef __MPI  
  ! Open like this only in sequential. Otherwize open with MPI-open
420
  IF ((etf_mem == 1) .AND. (ionode)) THEN
sponce's avatar
sponce committed
421 422 423 424 425
    ! open the .epmatwe file with the proper record length
    lrepmatw   = 2 * nbndsub * nbndsub * nrr_k * nmodes
    filint    = trim(prefix)//'.epmatwp'
    CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
  ENDIF
426
#endif
sponce's avatar
sponce committed
427 428
  ! 
  ! At this point, we will interpolate the Wannier rep to the Bloch rep 
sponce's avatar
sponce committed
429 430 431 432 433 434
  !
  IF ( epwread ) THEN
     !
     !  read all quantities in Wannier representation from file
     !  in parallel case all pools read the same file
     !
435
     CALL epw_read(nrr_k, nrr_q, nrr_g)
sponce's avatar
sponce committed
436 437
     !
  ELSE !if not epwread (i.e. need to calculate fmt file)
438
     ! 
439
     IF ((etf_mem == 1) .AND. (ionode)) THEN
440 441 442
       lrepmatw   = 2 * nbndsub * nbndsub * nrr_k * nmodes
       filint    = trim(prefix)//'.epmatwe'
       CALL diropn (iunepmatwe, 'epmatwe', lrepmatw, exst)
sponce's avatar
sponce committed
443
#ifdef __MPI       
444
       filint    = trim(prefix)//'.epmatwp'
445
       CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
sponce's avatar
sponce committed
446
#endif
447
     ENDIF
sponce's avatar
sponce committed
448
     !
449 450 451
     !xxq = 0.d0 
     !CALL loadumat &
     !     ( nbnd, nbndsub, nks, nkstot, xxq, cu, cuq, lwin, lwinq, exband )  
sponce's avatar
sponce committed
452 453 454 455 456
     !
     ! ------------------------------------------------------
     !   Bloch to Wannier transform
     ! ------------------------------------------------------
     !
Samuel Poncé's avatar
Samuel Poncé committed
457 458 459 460 461 462 463 464
     ALLOCATE( chw    ( nbndsub, nbndsub, nrr_k ),        &
               chw_ks ( nbndsub, nbndsub, nrr_k ),        &
               rdw    ( nmodes,  nmodes,  nrr_q ) )
     IF (vme) THEN 
       ALLOCATE( cvmew ( 3, nbndsub, nbndsub, nrr_k ) )
     ELSE
       ALLOCATE( cdmew ( 3, nbndsub, nbndsub, nrr_k ) )
     ENDIF
sponce's avatar
sponce committed
465 466
     ! 
     ! SP : Let the user chose. If false use files on disk
467
     IF (etf_mem == 0) THEN
sponce's avatar
sponce committed
468
       ALLOCATE(epmatwe ( nbndsub, nbndsub, nrr_k, nmodes, nqc))
469
       ALLOCATE (epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_g))
sponce's avatar
sponce committed
470 471
     ELSE
       ALLOCATE(epmatwe_mem ( nbndsub, nbndsub, nrr_k, nmodes))
Samuel Poncé's avatar
Samuel Poncé committed
472
       epmatwe_mem(:,:,:,:) = czero
sponce's avatar
sponce committed
473
     ENDIF
sponce's avatar
sponce committed
474 475 476 477
     !
     ! Hamiltonian
     !
     CALL hambloch2wan &
478
          ( nbnd, nbndsub, nks, nkstot, et, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw )
sponce's avatar
sponce committed
479 480 481 482
     !
     ! Kohn-Sham eigenvalues
     !
     IF (eig_read) THEN
483
       WRITE (stdout,'(5x,a)') "Interpolating MB and KS eigenvalues"
sponce's avatar
sponce committed
484
       CALL hambloch2wan &
485
            ( nbnd, nbndsub, nks, nkstot, et_ks, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw_ks )
sponce's avatar
sponce committed
486 487
     ENDIF
     !
Samuel Poncé's avatar
Samuel Poncé committed
488 489 490 491
     IF (vme) THEN 
       ! Transform of position matrix elements
       ! PRB 74 195118  (2006)
       CALL vmebloch2wan &
492
            ( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
Samuel Poncé's avatar
Samuel Poncé committed
493 494 495
     ELSE
       ! Dipole
       CALL dmebloch2wan &
496
            ( nbnd, nbndsub, nks, nkstot, dmec, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
Samuel Poncé's avatar
Samuel Poncé committed
497
     ENDIF
sponce's avatar
sponce committed
498
     !
Samuel Poncé's avatar
Samuel Poncé committed
499
     ! Dynamical Matrix
sponce's avatar
sponce committed
500
     !
501
     IF (.not. lifc) CALL dynbloch2wan &
502
                          ( nmodes, nqc, xqc, dynq, nrr_q, irvec_q, wslen_q )
sponce's avatar
sponce committed
503 504 505 506 507
     !
     !
     ! Electron-Phonon vertex (Bloch el and Bloch ph -> Wannier el and Bloch ph)
     !
     DO iq = 1, nqc
sponce's avatar
sponce committed
508 509 510 511 512
       !
       xxq = xqc (:, iq)
       !
       ! we need the cu again for the k+q points, we generate the map here
       !
513
       CALL loadumat ( nbnd, nbndsub, nks, nkstot, xxq, cu, cuq, lwin, lwinq, exband, w_centers )
sponce's avatar
sponce committed
514 515 516
       !
       DO imode = 1, nmodes
         !
517
         IF (etf_mem == 0) THEN 
sponce's avatar
sponce committed
518
           CALL ephbloch2wane &
519
             ( nbnd, nbndsub, nks, nkstot, xk, cu, cuq, &
520
             epmatq(:,:,:,imode,iq), nrr_k, irvec_k, wslen_k, epmatwe(:,:,:,imode,iq) )
sponce's avatar
sponce committed
521
         ELSE
sponce's avatar
sponce committed
522
           CALL ephbloch2wane &
523
             ( nbnd, nbndsub, nks, nkstot, xk, cu, cuq, &
524
             epmatq(:,:,:,imode,iq), nrr_k, irvec_k, wslen_k, epmatwe_mem(:,:,:,imode) )
sponce's avatar
sponce committed
525
           !
sponce's avatar
sponce committed
526 527 528 529
         ENDIF
         !
       ENDDO
       ! Only the master node writes 
530
       IF ((etf_mem == 1) .AND. (ionode)) THEN
sponce's avatar
sponce committed
531 532 533 534 535
         ! direct write of epmatwe for this iq 
         CALL rwepmatw ( epmatwe_mem, nbndsub, nrr_k, nmodes, iq, iunepmatwe, +1)       
         !   
       ENDIF   
       !
sponce's avatar
sponce committed
536 537 538 539
     ENDDO
     !
     ! Electron-Phonon vertex (Wannier el and Bloch ph -> Wannier el and Wannier ph)
     !
sponce's avatar
sponce committed
540 541
     ! Only master perform this task. Need to be parallelize in the future (SP)
     IF (ionode) THEN
542
       IF (etf_mem == 0) THEN
sponce's avatar
sponce committed
543
         CALL ephbloch2wanp &
544
           ( nbndsub, nmodes, xqc, nqc, irvec_k, irvec_g, nrr_k, nrr_g, epmatwe )
sponce's avatar
sponce committed
545 546
       ELSE
          CALL ephbloch2wanp_mem &
547
           ( nbndsub, nmodes, xqc, nqc, irvec_k, irvec_g, nrr_k, nrr_g, epmatwe_mem )
sponce's avatar
sponce committed
548 549 550 551
       ENDIF
     ENDIF
     !
     CALL mp_barrier(inter_pool_comm)
sponce's avatar
sponce committed
552 553
     !
     IF ( epwwrite ) THEN
554
        CALL epw_write(nrr_k, nrr_q, nrr_g, w_centers) 
555
        CALL epw_read(nrr_k, nrr_q, nrr_g) 
sponce's avatar
sponce committed
556 557 558 559 560
     ENDIF
     !
  ENDIF
  !
  IF ( ALLOCATED (epmatwe) ) DEALLOCATE (epmatwe)
sponce's avatar
sponce committed
561
  IF ( ALLOCATED (epmatwe_mem) ) DEALLOCATE (epmatwe_mem)
sponce's avatar
sponce committed
562 563 564 565 566
  IF ( ALLOCATED (epmatq) )  DEALLOCATE (epmatq)
  IF ( ALLOCATED (cu) )      DEALLOCATE (cu)
  IF ( ALLOCATED (cuq) )     DEALLOCATE (cuq)
  IF ( ALLOCATED (lwin) )    DEALLOCATE (lwin)
  IF ( ALLOCATED (lwinq) )   DEALLOCATE (lwinq)
Samuel Poncé's avatar
Samuel Poncé committed
567
  IF ( ALLOCATED (exband) )  DEALLOCATE (exband)
568 569 570 571 572
  IF (etf_mem == 1) THEN
    CLOSE(iunepmatwe, status = 'delete')
  ELSE
    CLOSE(iunepmatwe)
  ENDIF
573 574 575
#ifdef __MPI
  CLOSE(iunepmatwp)
#endif
sponce's avatar
sponce committed
576 577 578 579 580 581 582 583 584
  ! 
  ! Check Memory usage
  CALL system_mem_usage(valueRSS)
  ! 
  WRITE(stdout, '(a)' )             '     ==================================================================='
  WRITE(stdout, '(a,i10,a)' ) '     Memory usage:  VmHWM =',valueRSS(2)/1024,'Mb'
  WRITE(stdout, '(a,i10,a)' ) '                   VmPeak =',valueRSS(1)/1024,'Mb'
  WRITE(stdout, '(a)' )             '     ==================================================================='
  WRITE(stdout, '(a)' )             '     '
585
  
sponce's avatar
sponce committed
586
  !
587 588
  !  At this point, we will interpolate the Wannier rep to the Bloch rep 
  !  for electrons, phonons and the ep-matrix
sponce's avatar
sponce committed
589 590 591 592
  !
  !  need to add some sort of parallelization (on g-vectors?)  what
  !  else can be done when we don't ever see the wfcs??
  !
593 594
  CALL loadqmesh_serial
  CALL loadkmesh_para
sponce's avatar
sponce committed
595
  !
Samuel Poncé's avatar
Samuel Poncé committed
596 597 598 599 600 601 602 603 604 605
  ALLOCATE ( epmatwef( nbndsub, nbndsub, nrr_k, nmodes), &
             wf( nmodes,  nqf ),                         &
             etf( nbndsub, nkqf),                        &
             etf_ks( nbndsub, nkqf),                     &
             epmatf( nbndsub, nbndsub, nmodes),          &
             cufkk( nbndsub, nbndsub),                   &
             cufkq( nbndsub, nbndsub),                   & 
             uf( nmodes, nmodes),                        &
             bmatf( nbndsub, nbndsub),                   & 
             eps_rpa( nmodes) )
606 607
  !
  ! Need to be initialized
Samuel Poncé's avatar
Samuel Poncé committed
608
  etf_ks(:,:) = zero
609
  epmatf(:,:,:) = czero
sponce's avatar
sponce committed
610
  !
Samuel Poncé's avatar
Samuel Poncé committed
611 612 613 614 615 616 617
  ! allocate velocity and dipole matrix elements after getting grid size
  !
  IF (vme) THEN 
     ALLOCATE ( vmef(3, nbndsub, nbndsub, 2 * nkf) )
  ELSE
     ALLOCATE ( dmef(3, nbndsub, nbndsub, 2 * nkf) )
  ENDIF
sponce's avatar
sponce committed
618
  !
619 620
  ALLOCATE(cfac(nrr_k,dims,dims))
  ALLOCATE(cfacq(nrr_k,dims,dims))
621
  ALLOCATE(rdotk(nrr_k))
622
  ALLOCATE(rdotk2(nrr_k))
623 624
  ! This is simply because dgemv take only real number (not integer)
  ALLOCATE(irvec_r(3,nrr_k))
625
  irvec_r = REAL(irvec_k,KIND=dp)
626
  ! 
sponce's avatar
sponce committed
627 628 629 630 631 632 633 634 635 636 637 638 639
  ! ------------------------------------------------------
  ! Hamiltonian : Wannier -> Bloch (preliminary)
  ! ------------------------------------------------------
  !
  ! We here perform a preliminary interpolation of the hamiltonian
  ! in order to determine the fermi window ibndmin:ibndmax for later use.
  ! We will interpolate again afterwards, for each k and k+q separately
  !
  xxq = 0.d0
  !
  ! nkqf is the number of kpoints in the pool
  !
  DO ik = 1, nkqf
640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668
    !
    xxk = xkf (:, ik)
    !
    IF ( 2*(ik/2).eq.ik ) THEN
      !
      !  this is a k+q point : redefine as xkf (:, ik-1) + xxq
      !
      CALL cryst_to_cart ( 1, xxq, at,-1 )
      xxk = xkf (:, ik-1) + xxq
      CALL cryst_to_cart ( 1, xxq, bg, 1 )
      !
    ENDIF
    !
    ! SP: Compute the cfac only once here since the same are use in both hamwan2bloch and dmewan2bloch
    ! + optimize the 2\pi r\cdot k with Blas
    CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xxk, 1, 0.0_DP, rdotk, 1 )
    ! 
    DO iw=1, dims
      DO iw2=1, dims
        DO ir=1, nrr_k
          IF (ndegen_k(ir,iw2,iw) > 0 ) &
            cfac(ir,iw2,iw) = exp( ci*rdotk(ir) ) / ndegen_k(ir,iw2,iw)
        ENDDO
      ENDDO
    ENDDO
    ! 
    CALL hamwan2bloch &
         ( nbndsub, nrr_k, cufkk, etf(:, ik), chw, cfac, dims)
    !
sponce's avatar
sponce committed
669 670 671 672 673 674 675 676 677 678 679
  ENDDO
  !
  ! 27/06/2012 RM
  ! in the case when a random or uniform fine k-mesh is used
  ! calculate the Fermi level corresponding to the fine k-mesh 
  ! this Fermi level is then used as a reference in fermiwindow 
  ! 06/05/2014 CV
  ! calculate the Fermi level corresponding to the fine k-mesh
  ! or read it from input (Fermi level from the coarse grid 
  ! may be wrong or inaccurate)
  !
680
  WRITE(stdout,'(/5x,a,f10.6,a)') 'Fermi energy coarse grid = ', ef * ryd2ev, ' eV'
sponce's avatar
sponce committed
681 682 683 684
  !
  IF( efermi_read ) THEN
     !
     ef = fermi_energy
685 686
     WRITE(stdout,'(/5x,a)') repeat('=',67)
     WRITE(stdout, '(/5x,a,f10.6,a)') &
sponce's avatar
sponce committed
687
         'Fermi energy is read from the input file: Ef = ', ef * ryd2ev, ' eV'
688
     WRITE(stdout,'(/5x,a)') repeat('=',67)
689
     !
690 691 692 693 694 695 696 697 698 699
     ! SP: even when reading from input the number of electron needs to be correct
     already_skipped = .false.
     IF ( nbndskip .gt. 0 ) THEN
        IF ( .not. already_skipped ) THEN
           IF ( noncolin ) THEN
              nelec = nelec - one * nbndskip
           ELSE
              nelec = nelec - two * nbndskip
           ENDIF
           already_skipped = .true.
700 701
           WRITE(stdout,'(/5x,"Skipping the first ",i4," bands:")') nbndskip
           WRITE(stdout,'(/5x,"The Fermi level will be determined with ",f9.5," electrons")') nelec
702
        ENDIF
703 704
     ENDIF
     !      
sponce's avatar
sponce committed
705 706
  ELSEIF( band_plot ) THEN 
     !
707
     WRITE(stdout,'(/5x,a)') repeat('=',67)
sponce's avatar
sponce committed
708
     WRITE(stdout, '(/5x,"Fermi energy corresponds to the coarse k-mesh")')
709
     WRITE(stdout,'(/5x,a)') repeat('=',67) 
sponce's avatar
sponce committed
710 711 712 713 714 715 716 717 718 719 720 721 722 723
     !
  ELSE 
     ! here we take into account that we may skip bands when we wannierize
     ! (spin-unpolarized)
     ! RM - add the noncolin case
     already_skipped = .false.
     IF ( nbndskip .gt. 0 ) THEN
        IF ( .not. already_skipped ) THEN
           IF ( noncolin ) THEN 
              nelec = nelec - one * nbndskip
           ELSE
              nelec = nelec - two * nbndskip
           ENDIF
           already_skipped = .true.
724 725
           WRITE(stdout,'(/5x,"Skipping the first ",i4," bands:")') nbndskip
           WRITE(stdout,'(/5x,"The Fermi level will be determined with ",f9.5," electrons")') nelec
sponce's avatar
sponce committed
726 727 728 729 730 731 732
        ENDIF
     ENDIF
     !
     ! Fermi energy
     !  
     ! since wkf(:,ikq) = 0 these bands do not bring any contribution to Fermi level
     !  
733
     efnew = efermig(etf, nbndsub, nkqf, nelec, wkf, degaussw, ngaussw, 0, isk)
sponce's avatar
sponce committed
734
     !
735
     WRITE(stdout, '(/5x,a,f10.6,a)') &
sponce's avatar
sponce committed
736 737 738 739 740
         'Fermi energy is calculated from the fine k-mesh: Ef = ', efnew * ryd2ev, ' eV'
     !
     ! if 'fine' Fermi level differs by more than 250 meV, there is probably something wrong
     ! with the wannier functions, or 'coarse' Fermi level is inaccurate
     IF (abs(efnew - ef) * ryd2eV .gt. 0.250d0 .and. (.not.eig_read) ) &
741 742
        WRITE(stdout,'(/5x,a)') 'Warning: check if difference with Fermi level fine grid makes sense'
     WRITE(stdout,'(/5x,a)') repeat('=',67)
sponce's avatar
sponce committed
743
     !
Samuel Poncé's avatar
Samuel Poncé committed
744
     ef = efnew
sponce's avatar
sponce committed
745 746 747 748 749 750 751 752
     !
  ENDIF
  !
  ! identify the bands within fsthick from the Fermi level
  ! (in shuffle mode this actually does not depend on q)
  !
  !
  CALL fermiwindow
753 754 755 756 757 758 759 760 761 762
  ! 
  ! Define it only once for the full run. 
  CALL fkbounds( nkqtotf/2, lower_bnd, upper_bnd )
  ! 
  ! Re-order the k-point according to weather they are in or out of the fshick
  ! windows
  IF (iterative_bte .and. mp_mesh_k) THEN
    CALL load_rebal() 
  ENDIF
  !
sponce's avatar
sponce committed
763 764 765 766 767
  !  xqf must be in crystal coordinates
  !
  ! this loops over the fine mesh of q points.
  ! ---------------------------------------------------------------------------------------
  ! ---------------------------------------------------------------------------------------
768 769 770 771 772 773 774 775 776 777
  IF (lifc) THEN
    !
    ! build the WS cell corresponding to the force constant grid
    !
    atws(:,1) = at(:,1)*DBLE(nq1)
    atws(:,2) = at(:,2)*DBLE(nq2)
    atws(:,3) = at(:,3)*DBLE(nq3)
    ! initialize WS r-vectors
    CALL wsinit(rws,nrwsx,nrws,atws)
  ENDIF
778 779 780 781 782 783 784 785 786 787 788
  !
  ! Open the ephmatwp file here
#if defined(__MPI)
  IF (etf_mem == 1) then
    ! Check for directory given by "outdir"
    !      
    filint = trim(tmp_dir)//trim(prefix)//'.epmatwp1'
    CALL MPI_FILE_OPEN(world_comm,filint,MPI_MODE_RDONLY,MPI_INFO_NULL,iunepmatwp2,ierr)
    IF( ierr /= 0 ) CALL errore( 'ephwann_shuffle', 'error in MPI_FILE_OPEN',1 )
  ENDIF
#endif
789
  !
790 791 792 793 794 795 796
  ! get the size of the matrix elements stored in each pool
  ! for informational purposes.  Not necessary
  !
  CALL mem_size(ibndmin, ibndmax, nmodes, nkf)
  !
  ! Fine mesh set of g-matrices.  It is large for memory storage
  ALLOCATE ( epf17 (ibndmax-ibndmin+1, ibndmax-ibndmin+1, nmodes, nkf) )
797
  ALLOCATE ( etf_all ( ibndmax-ibndmin+1, nkqtotf/2 ) )
798 799
  ALLOCATE ( inv_tau_all (nstemp, ibndmax-ibndmin+1, nkqtotf/2) )
  ALLOCATE ( zi_allvb (nstemp, ibndmax-ibndmin+1, nkqtotf/2) )
800 801 802 803
  epf17(:,:,:,:)     = czero 
  etf_all(:,:)       = zero
  inv_tau_all(:,:,:) = zero
  zi_allvb(:,:,:)    = zero
804 805 806 807
  ! 
  IF (int_mob .AND. carrier) THEN
    ALLOCATE ( inv_tau_allcb (nstemp, ibndmax-ibndmin+1, nkqtotf/2) )
    ALLOCATE ( zi_allcb (nstemp, ibndmax-ibndmin+1, nkqtotf/2) )
808 809
    inv_tau_allcb(:,:,:) = zero
    zi_allcb(:,:,:)      = zero
810
  ENDIF
811 812 813 814 815 816 817
  ! 
  ! ------------------------------------------------
  ! The IBTE implement works in two steps
  ! 1) compute the dominant scattering rates and store them to file
  ! 2) read them from file and solve the IBTE where all important element are in memory
  ! ------------------------------------------------
  !  
818
  ! Initialization and restart when doing IBTE
819
  IF (iterative_bte) THEN
820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837
998 continue ! Continue after all scattering rates have been computed in print_ibte
    IF (epmatkqread) THEN
      ! 
      ALLOCATE( vkk_all( 3, ibndmax-ibndmin+1, nkqtotf/2 ) )
      ALLOCATE( wkf_all( nkqtotf/2 ) )
      ! 
      CALL iter_restart(etf_all, wkf_all, vkk_all, ind_tot, ind_totcb, ef0, efcb)
      ! 
      DEALLOCATE(vkk_all)
      DEALLOCATE(wkf_all)
      DEALLOCATE(etf_all)
      GOTO 999
      ! 
    ELSE ! epmatkqread
      !  
      ! Open the required files
      CALL iter_open(ind_tot, ind_totcb, lrepmatw2, lrepmatw4, lrepmatw5, lrepmatw6) 
      ! 
Samuel Poncé's avatar
Samuel Poncé committed
838
    ENDIF
839
  ENDIF 
840
  ! 
841 842 843 844 845
  ! -----------------------------------------------------------------------
  ! Determines which q-points falls within the fsthick windows
  ! Store the result in the selecq.fmt file 
  ! If the file exists, automatically restart from the file
  ! -----------------------------------------------------------------------
846 847 848 849 850 851 852 853 854 855 856
  ! 
  ! Check if the grids are homogeneous and commensurate
  homogeneous = .FALSE.
  IF ( (nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0) .AND. &
       (nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0) .AND. &
       (MOD(nkf1,nqf1) == 0) .AND. (MOD(nkf2,nqf2) == 0) .AND. (MOD(nkf3,nqf3) == 0) ) THEN
    homogeneous = .TRUE.
  ELSE
    homogeneous = .FALSE.
  ENDIF
  ! 
857
  totq = 0
858 859 860 861 862
  ! Check if the file has been pre-computed
  IF (mpime == ionode_id) THEN
    INQUIRE(FILE='selecq.fmt',EXIST=exst)
  ENDIF
  CALL mp_bcast(exst, ionode_id, world_comm)
863
  ! 
864 865 866 867 868 869 870 871 872
  IF (exst) THEN
    IF (selecqread) THEN
      WRITE(stdout,'(5x,a)')' '
      WRITE(stdout,'(5x,a)')'Reading selecq.fmt file. '
      CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
    ELSE 
      WRITE(stdout,'(5x,a)')' '
      WRITE(stdout,'(5x,a)')'A selecq.fmt file was found but re-created because selecqread == .false. '
      CALL qwindow(.FALSE., nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
873
    ENDIF
874 875 876 877 878
  ELSE ! exst
    IF (selecqread) THEN
      CALL errore( 'ephwann_shuffle', 'Variable selecqread == .true. but file selecq.fmt not found.',1 ) 
    ELSE
      CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
879
    ENDIF
880 881 882 883
  ENDIF
  ! 
  WRITE(stdout,'(5x,a,i8,a)')'We only need to compute ',totq, ' q-points'
  WRITE(stdout,'(5x,a)')' '
884 885 886 887
  ! 
  ! -----------------------------------------------------------------------
  ! Possible restart during step 1) 
  ! -----------------------------------------------------------------------
888 889 890
  iq_restart = 1
  first_cycle = .FALSE.
  first_time = .TRUE.
891 892
  ! 
  ! Restart in SERTA case or self-energy case
893 894
  IF (restart) THEN
    IF ( elecselfen ) THEN
895 896 897
      ALLOCATE( sigmar_all(ibndmax-ibndmin+1, nkqtotf/2) )
      ALLOCATE( sigmai_all(ibndmax-ibndmin+1, nkqtotf/2) )
      ALLOCATE( zi_all(ibndmax-ibndmin+1, nkqtotf/2) )
898 899
      sigmar_all(:,:) = zero
      sigmai_all(:,:) = zero
900 901
      zi_all(:,:)     = zero
      CALL electron_read(iq_restart, totq, nkqtotf/2, sigmar_all, sigmai_all, zi_all)
902
    ENDIF
903 904 905
    IF ( scattering ) THEN
      IF (int_mob .AND. carrier) THEN
        ! Here inv_tau_all and inv_tau_allcb gets updated
906
        CALL tau_read(iq_restart, totq, nkqtotf/2, .TRUE.)
907 908
      ELSE
        ! Here inv_tau_all gets updated
909
        CALL tau_read(iq_restart, totq, nkqtotf/2, .FALSE.)
910 911
      ENDIF
    ENDIF
912 913 914 915 916
    !
    ! If you restart from reading a file. This prevent 
    ! the case were you restart but the file does not exist
    IF (iq_restart > 1) first_cycle = .TRUE.
    ! 
917
  ENDIF ! restart
918 919 920
  ! 
  ! Scatread assumes that you alread have done the full q-integration
  ! We just do one loop to get interpolated eigenenergies.  
921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940
  IF(scatread) iq_restart = totq -1
  ! 
  ! Restart in IBTE case
  IF (iterative_bte) THEN
    IF (mpime == ionode_id) THEN
      INQUIRE(FILE='restart_ibte.fmt',EXIST=exst)
    ENDIF
    CALL mp_bcast(exst, ionode_id, world_comm)
    ! 
    IF (exst) THEN
      IF (mpime.eq.ionode_id) THEN
        OPEN(unit=iunrestart,file='restart_ibte.fmt',status='old',iostat=ios)
        READ (iunrestart,*) iq_restart
        READ (iunrestart,*) ind_tot
        READ (iunrestart,*) ind_totcb
        READ (iunrestart,*) lrepmatw2
        READ (iunrestart,*) lrepmatw4
        READ (iunrestart,*) lrepmatw5
        READ (iunrestart,*) lrepmatw6
        CLOSE(iunrestart)
941
      ENDIF
942
      CALL mp_bcast(iq_restart, ionode_id, world_comm )
943
#if defined(__MPI)
944 945 946 947 948 949
      CALL MPI_BCAST( ind_tot,   1, MPI_OFFSET, ionode_id, world_comm, ierr)
      CALL MPI_BCAST( ind_totcb, 1, MPI_OFFSET, ionode_id, world_comm, ierr)
      CALL MPI_BCAST( lrepmatw2, 1, MPI_OFFSET, ionode_id, world_comm, ierr)
      CALL MPI_BCAST( lrepmatw4, 1, MPI_OFFSET, ionode_id, world_comm, ierr)
      CALL MPI_BCAST( lrepmatw5, 1, MPI_OFFSET, ionode_id, world_comm, ierr)
      CALL MPI_BCAST( lrepmatw6, 1, MPI_OFFSET, ionode_id, world_comm, ierr)
950 951 952 953 954 955 956 957
#else
      CALL mp_bcast( ind_tot,   ionode_id, world_comm )
      CALL mp_bcast( ind_totcb, ionode_id, world_comm )
      CALL mp_bcast( lrepmatw2, ionode_id, world_comm )
      CALL mp_bcast( lrepmatw4, ionode_id, world_comm )
      CALL mp_bcast( lrepmatw5, ionode_id, world_comm )
      CALL mp_bcast( lrepmatw6, ionode_id, world_comm )
#endif
958
      IF( ierr /= 0 ) CALL errore( 'ephwann_shuffle', 'error in MPI_BCAST',1 )
959
      ! 
960 961 962 963 964
      ! Now, the iq_restart point has been done, so we need to do the next one except if last
      !IF (iq_restart /= totq) iq_restart = iq_restart + 1
      ! Now, the iq_restart point has been done, so we need to do the next 
      iq_restart = iq_restart + 1
      WRITE(stdout,'(5x,a,i8,a)')'We restart from ',iq_restart, ' q-points'
965
      ! 
966 967 968 969 970 971
    ENDIF ! exst
  ENDIF
  ! -----------------------------------------------------------------------------
  ! 
  DO iqq = iq_restart, totq
    ! This needs to be uncommented. 
972
    epf17(:,:,:,:) = czero
973 974 975 976 977 978 979 980 981
    ! 
    iq = selecq(iqq)
    !   
    CALL start_clock ( 'ep-interp' )
    !
    ! In case of big calculation, show progression of iq (especially usefull when
    ! elecselfen = true as nothing happen during the calculation otherwise. 
    !
    IF ( .not. phonselfen) THEN 
982 983
      IF (MOD(iqq,100) == 0) THEN
        WRITE(stdout, '(5x,a,i10,a,i10)' ) 'Progression iq (fine) = ',iqq,'/',totq
984
      ENDIF
985 986
    ENDIF
    !
987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
    xxq = xqf (:, iq)
    !
    ! ------------------------------------------------------
    ! dynamical matrix : Wannier -> Bloch
    ! ------------------------------------------------------
    !
    IF (.not. lifc) THEN
      CALL dynwan2bloch &
          ( nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2 )
    ELSE
      CALL dynifc2blochf ( nmodes, rws, nrws, xxq, uf, w2 )
    ENDIF
    !
    ! ...then take into account the mass factors and square-root the frequencies...
    !
    DO nu = 1, nmodes
      !
      ! wf are the interpolated eigenfrequencies
      ! (omega on fine grid)
      !
      IF ( w2 (nu) .gt. 0.d0 ) THEN
        wf(nu,iq) =  sqrt(abs( w2 (nu) ))
      ELSE
        wf(nu,iq) = -sqrt(abs( w2 (nu) ))
      ENDIF
      !
      DO mu = 1, nmodes
        na = (mu - 1) / 3 + 1
        uf (mu, nu) = uf (mu, nu) / sqrt(amass(ityp(na)))
      ENDDO
    ENDDO
    !
    ! --------------------------------------------------------------
    ! epmat : Wannier el and Wannier ph -> Wannier el and Bloch ph
    ! --------------------------------------------------------------
    !
    !DBSP
    !CALL start_clock ( 'cl2' )
    IF (.NOT. longrange) THEN
      CALL ephwan2blochp &
          ( nmodes, xxq, irvec_g, ndegen_g, nrr_g, uf, epmatwef, nbndsub, nrr_k, dims, dims2 )
    ENDIF
    !CALL stop_clock ( 'cl2' )
    !
    !
    !  number of k points with a band on the Fermi surface
    fermicount = 0
    !
    IF (lscreen) THEN
      IF (scr_typ == 0) CALL rpa_epsilon (xxq, wf(:,iq), nmodes, epsi, eps_rpa)
      IF (scr_typ == 1) CALL tf_epsilon (xxq, nmodes, epsi, eps_rpa)
    ENDIF
    ! 
    ! this is a loop over k blocks in the pool
    ! (size of the local k-set)
    DO ik = 1, nkf
      !
      ! xkf is assumed to be in crys coord
      !
      ikk = 2 * ik - 1
      ikq = ikk + 1
      !
      xkk = xkf(:, ikk)
      xkq = xkk + xxq
      !
      CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
      CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkq, 1, 0.0_DP, rdotk2, 1 )
      !
      IF (use_ws) THEN
        DO iw=1, dims
          DO iw2=1, dims
            DO ir = 1, nrr_k
              IF (ndegen_k(ir,iw2,iw) > 0) THEN
                cfac(ir,iw2,iw)  = exp( ci*rdotk(ir) ) / ndegen_k(ir,iw2,iw)
                cfacq(ir,iw2,iw) = exp( ci*rdotk2(ir) ) / ndegen_k(ir,iw2,iw)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ELSE 
        cfac(:,1,1)   = exp( ci*rdotk(:) ) / ndegen_k(:,1,1)
        cfacq(:,1,1)  = exp( ci*rdotk2(:) ) / ndegen_k(:,1,1)
      ENDIF
      !
      ! ------------------------------------------------------        
      ! hamiltonian : Wannier -> Bloch 
      ! ------------------------------------------------------
      !
      ! Kohn-Sham first, then get the rotation matricies for following interp.
      IF (eig_read) THEN
         CALL hamwan2bloch &
           ( nbndsub, nrr_k, cufkk, etf_ks(:, ikk), chw_ks, cfac, dims)
         CALL hamwan2bloch &
           ( nbndsub, nrr_k, cufkq, etf_ks(:, ikq), chw_ks, cfacq, dims)
      ENDIF
      !
      CALL hamwan2bloch &
           ( nbndsub, nrr_k, cufkk, etf(:, ikk), chw, cfac, dims)
      CALL hamwan2bloch &
           ( nbndsub, nrr_k, cufkq, etf(:, ikq), chw, cfacq, dims)
      !
      IF (vme) THEN
1089
         !
1090 1091 1092
         ! ------------------------------------------------------
         !  velocity: Wannier -> Bloch
         ! ------------------------------------------------------
1093
         !
1094 1095 1096 1097 1098
         IF (eig_read) THEN
            CALL vmewan2bloch &
                 ( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw_ks, cfac, dims )
            CALL vmewan2bloch &
                 ( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw_ks, cfacq, dims )
1099
         ELSE
1100 1101 1102 1103
            CALL vmewan2bloch &
                 ( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw, cfac, dims )
            CALL vmewan2bloch &
                 ( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw, cfacq, dims )
1104
         ENDIF
1105
      ELSE
1106
         !
1107 1108
         ! ------------------------------------------------------
         !  dipole: Wannier -> Bloch
1109
         ! ------------------------------------------------------
1110
         !
1111 1112 1113 1114
         CALL dmewan2bloch &
              ( nbndsub, nrr_k, cufkk, dmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), cfac, dims)
         CALL dmewan2bloch &
              ( nbndsub, nrr_k, cufkq, dmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), cfacq, dims)
1115
         !
1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139
      ENDIF
      !
      IF (.NOT. scatread) THEN
        ! interpolate only when (k,k+q) both have at least one band 
        ! within a Fermi shell of size fsthick 
        !
        IF ( (( minval ( abs(etf(:, ikk) - ef) ) < fsthick ) .and. & 
              ( minval ( abs(etf(:, ikq) - ef) ) < fsthick )) ) THEN
          !
          !  fermicount = fermicount + 1
          !
          ! --------------------------------------------------------------
          ! epmat : Wannier el and Bloch ph -> Bloch el and Bloch ph
          ! --------------------------------------------------------------
          !
          ! SP: Note: In case of polar materials, computing the long-range and short-range term 
          !     separately might help speed up the convergence. Indeed the long-range term should be 
          !     much faster to compute. Note however that the short-range term still contains a linear
          !     long-range part and therefore could still be a bit more difficult to converge than 
          !     non-polar materials. 
          ! 
          IF (longrange) THEN
            !      
            epmatf(:,:,:) = czero
1140
            !
1141
          ELSE
1142
            !
1143 1144 1145
            epmatf(:,:,:) = czero
            CALL ephwan2bloch &
              ( nbndsub, nrr_k, epmatwef, cufkk, cufkq, epmatf, nmodes, cfac, dims )
1146
            !
1147 1148 1149
          ENDIF
          !
          IF (lpolar) THEN
1150
            !
1151 1152
            CALL compute_umn_f( nbndsub, cufkk, cufkq, bmatf )
            !
1153
            IF ( (abs(xxq(1)) > eps8) .or. (abs(xxq(2)) > eps8) .or. (abs(xxq(3)) > eps8) ) THEN
1154 1155 1156 1157 1158 1159 1160
              !      
              CALL cryst_to_cart (1, xxq, bg, 1)
              CALL rgd_blk_epw_fine(nq1, nq2, nq3, xxq, uf, epmatf, &
                                    nmodes, epsi, zstar, bmatf, one)
              CALL cryst_to_cart (1, xxq, at, -1)
              !
            ENDIF
1161
            !
1162
          ENDIF
1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181
          ! 
          ! Store epmatf in memory
          !
          DO jbnd = ibndmin, ibndmax
            DO ibnd = ibndmin, ibndmax
              ! 
              IF (lscreen) THEN
                 epf17(ibnd-ibndmin+1,jbnd-ibndmin+1,:,ik) = epmatf(ibnd,jbnd,:) / eps_rpa(:)
              ELSE
                 epf17(ibnd-ibndmin+1,jbnd-ibndmin+1,:,ik) = epmatf(ibnd,jbnd,:)
              ENDIF
              !
            ENDDO
          ENDDO
          !
        ENDIF
      ENDIF ! scatread 
      !
    ENDDO  ! end loop over k points
1182
    !
1183 1184 1185 1186 1187 1188 1189 1190 1191
    IF (prtgkk     ) CALL print_gkk( iq )
    IF (phonselfen ) CALL selfen_phon_q( iqq, iq, totq )
    IF (elecselfen ) CALL selfen_elec_q( iqq, iq, totq, first_cycle )
    IF (plselfen .and. .not. vme ) CALL selfen_pl_q( iqq, iq, totq )
    IF (nest_fn    ) CALL nesting_fn_q( iqq, iq )
    IF (specfun_el ) CALL spectral_func_q( iqq, iq, totq )
    IF (specfun_ph ) CALL spectral_func_ph( iqq, iq, totq )
    IF (specfun_pl .and. .not. vme ) CALL spectral_func_pl_q( iqq, iq, totq )
    IF (ephwrite) THEN
1192
      IF ( iq == 1 ) THEN 
1193 1194 1195
         CALL kmesh_fine
         CALL kqmap_fine
      ENDIF
1196 1197
      CALL write_ephmat( iq ) 
      CALL count_kpoints( iq )
1198
    ENDIF
sponce's avatar
sponce committed
1199
    ! 
1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223
    IF (.NOT. scatread) THEN
      ! 
      ! Indirect absorption ---------------------------------------------------------
      ! If Indirect absortpion, keep unshifted values:
      IF ( lindabs .AND. .NOT. scattering ) etf_ks(:,:) = etf(:,:)
      ! 
      ! Apply a scissor shift to CBM if required by user
      ! The shift is apply to k and k+q
      IF (ABS(scissor) > eps6) THEN
        IF ( noncolin ) THEN
          icbm = FLOOR(nelec/1.0d0) +1
        ELSE
          icbm = FLOOR(nelec/2.0d0) +1
        ENDIF
        !
        DO ik = 1, nkf
          ikk = 2 * ik - 1
          ikq = ikk + 1
          DO ibnd = icbm, nbndsub
            ! 
            etf (ibnd, ikk) = etf (ibnd, ikk) + scissor
            etf (ibnd, ikq) = etf (ibnd, ikq) + scissor
          ENDDO
        ENDDO
1224
        IF ( iq == 1 ) THEN
1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333
          WRITE(stdout, '(5x,"Applying a scissor shift of ",f9.5," eV to the conduction states")' ) scissor * ryd2ev
        ENDIF
      ENDIF
      !  
      ! Indirect absorption
      IF ( lindabs .AND. .NOT. scattering )  CALL indabs(iq)  
      ! 
      ! Conductivity ---------------------------------------------------------
      IF (scattering) THEN
        !   
        ! If we want to compute intrinsic mobilities, call fermicarrier to 
        ! correctly positionned the ef0 level.
        ! This is only done once for iq = 0 
        IF ( iqq == iq_restart ) THEN
          ! 
          DO itemp = 1, nstemp
            ! 
            etemp = transp_temp(itemp) 
            WRITE(stdout, '(/5x,"Temperature ",f8.3," K")' ) etemp * ryd2ev / kelvin2eV
            ! 
            ! Small gap semiconductor. Computes intrinsic mobility by placing 
            ! the Fermi level such that carrier density is equal for electron and holes
            IF (int_mob .AND. .NOT. carrier) THEN               
              !
              ef0(itemp) = fermicarrier( etemp )
              WRITE(stdout, '(5x,"Mobility Fermi level ",f10.6," eV")' )  ef0(itemp) * ryd2ev  
              ! We only compute 1 Fermi level so we do not need the other
              efcb(itemp) = 0
              !   
            ENDIF
            ! 
            ! Large bandgap semiconductor. Place the gap at the value ncarrier.
            ! The user want both VB and CB mobilities. 
            IF (int_mob .AND. carrier) THEN
              ! 
              ncarrier = - ABS(ncarrier) 
              ef0(itemp) = fermicarrier( etemp )
              WRITE(stdout, '(5x,"Mobility VB Fermi level ",f10.6," eV")' )  ef0(itemp) * ryd2ev 
              ! 
              ncarrier = ABS(ncarrier) 
              efcb(itemp) = fermicarrier( etemp )
              WRITE(stdout, '(5x,"Mobility CB Fermi level ",f10.6," eV")' )  efcb(itemp) * ryd2ev
              !  
            ENDIF   
            ! 
            ! User decide the carrier concentration and choose to only look at VB or CB  
            IF (.NOT. int_mob .AND. carrier) THEN
              ! SP: Determination of the Fermi level for intrinsic or doped carrier 
              !     One also need to apply scissor before calling it.
              ! 
              ef0(itemp) = fermicarrier( etemp )               
              WRITE(stdout, '(5x,"Mobility Fermi level ",f10.6," eV")' )  ef0(itemp) * ryd2ev
              ! We only compute 1 Fermi level so we do not need the other
              efcb(itemp) = 0
              ! 
            ENDIF
            ! 
            IF (.NOT. int_mob .AND. .NOT. carrier ) THEN
              IF ( efermi_read ) THEN
                !
                ef0(itemp) = fermi_energy
                !
              ELSE !SP: This is added for efficiency reason because the efermig routine is slow
                ef0(itemp) = efnew
              ENDIF
              ! We only compute 1 Fermi level so we do not need the other
              efcb(itemp) = 0
              !  
            ENDIF
            ! 
          ENDDO
          !
          ! 
        ENDIF ! iqq=0
        !   
        IF ( .NOT. iterative_bte ) THEN
          CALL scattering_rate_q( iqq, iq, totq, ef0, efcb, first_cycle )
          ! Computes the SERTA mobility
          !IF (iq == nqf) CALL transport_coeffs (ef0,efcb)
          IF (iqq == totq) CALL transport_coeffs (ef0,efcb)
        ENDIF
        ! 
        IF (iterative_bte) THEN
          CALL print_ibte( iqq, iq, totq, ef0, efcb, first_cycle, ind_tot, ind_totcb, lrepmatw2,&
                           lrepmatw4, lrepmatw5, lrepmatw6 )
          !  
          ! Finished, now compute SERTA and IBTE mobilities
          IF (iqq == totq) THEN
            WRITE(stdout, '(5x,a)')' '
            WRITE(stdout, '(5x,"epmatkqread automatically changed to .true. as all scattering have been computed.")')
            WRITE(stdout, '(5x,a)')' '
            ! close files
            CALL iter_close()
            !   
            epmatkqread = .true.
            GOTO 998
          ENDIF  
        ENDIF
        ! 
      ENDIF ! scattering
      ! --------------------------------------       
      !
      CALL stop_clock ( 'ep-interp' )
      !
    ENDIF ! scatread
  ENDDO  ! end loop over q points
  !
  ! close files
  CALL iter_close() 
sponce's avatar
sponce committed
1334
  ! 
1335 1336 1337 1338 1339 1340 1341 1342 1343
  ! Check Memory usage
  CALL system_mem_usage(valueRSS)
  ! 
  WRITE(stdout, '(a)' )             '     ==================================================================='
  WRITE(stdout, '(a,i10,a)' ) '     Memory usage:  VmHWM =',valueRSS(2)/1024,'Mb'
  WRITE(stdout, '(a,i10,a)' ) '                   VmPeak =',valueRSS(1)/1024,'Mb'
  WRITE(stdout, '(a)' )             '     ==================================================================='
  WRITE(stdout, '(a)' )
  ! 
sponce's avatar
sponce committed
1344 1345 1346 1347 1348 1349 1350 1351
  ! ---------------------------------------------------------------------------------------
  ! ---------------------------------------------------------------------------------------  
  !
  ! SP: Added lambda and phonon lifetime writing to file.
  ! 
  CALL mp_barrier(inter_pool_comm)
  IF (mpime.eq.ionode_id) THEN
    !
1352
    IF (phonselfen) THEN
sponce's avatar
sponce committed
1353 1354 1355
      OPEN(unit=lambda_phself,file='lambda.phself')
      WRITE(lambda_phself, '(/2x,a/)') '#Lambda phonon self-energy'
      WRITE(lambda_phself, *) '#Modes     ',(imode, imode=1,nmodes)
1356
      DO iqq = 1, nqtotf
sponce's avatar
sponce committed
1357
          !
1358 1359
          !myfmt = "(*(3x,E15.5))"  This does not work with PGI
        myfmt = "(1000(3x,E15.5))"
1360 1361
        WRITE(lambda_phself,'(i9,4x)',advance='no') iqq
        WRITE(lambda_phself, fmt=myfmt) (REAL(lambda_all(imode,iqq,1)),imode=1,nmodes)
sponce's avatar
sponce committed
1362 1363 1364 1365
          !
      ENDDO
      CLOSE(lambda_phself)
      OPEN(unit=linewidth_phself,file='linewidth.phself')
1366 1367
      WRITE(linewidth_phself, '(a)') '# Phonon frequency and phonon lifetime in meV '
      WRITE(linewidth_phself,'(a)') '# Q-point  Mode   Phonon freq (meV)   Phonon linewidth (meV)'
1368
      DO iqq = 1, nqtotf
sponce's avatar
sponce committed
1369
        !
1370
        DO imode=1, nmodes
Samuel Poncé's avatar