*** /home/giannozz/espresso-old/espresso-5.4.0/PW/src/pwcom.f90	2016-04-24 23:19:18.000000000 +0200
--- PW/src/pwcom.f90	2016-06-21 14:04:55.342280534 +0200
***************
*** 7,12 ****
--- 7,13 ----
  !
  !--------------------------------------------------------------------------
  !
+ !
  MODULE klist
    !
    ! ... The variables for the k-points
***************
*** 14,19 ****
--- 15,21 ----
    USE kinds,      ONLY : DP
    USE parameters, ONLY : npk
    !
+   IMPLICIT NONE
    SAVE
    !
    CHARACTER (len=32) :: &
***************
*** 44,52 ****
         two_fermi_energies ! if .TRUE.: nelup and neldw set ef_up and ef_dw
                            ! separately
    !
  END MODULE klist
  !
- !
  MODULE lsda_mod
    !
    ! ... The variables needed for the lsda calculation
--- 46,88 ----
         two_fermi_energies ! if .TRUE.: nelup and neldw set ef_up and ef_dw
                            ! separately
    !
+ CONTAINS
+   !
+   SUBROUTINE init_igk ( npwx, ngm, g, gcutw )
+     !
+     ! ... Initialize indices igk_k and number of plane waves per k-point:
+     ! ...    (k_ik+G)_i = k_ik+G_igk,   i=1,ngk(ik), igk=igk_k(i,ik)
+     !
+     INTEGER, INTENT (IN) :: npwx, ngm
+     REAL(dp), INTENT(IN) :: gcutw, g(3,ngm)
+     !
+     REAL(dp), ALLOCATABLE :: gk (:)
+     INTEGER :: ik
+     !
+ 
+     IF(.NOT.ALLOCATED(igk_k)) ALLOCATE ( igk_k(npwx,nks))
+     IF(.NOT.ALLOCATED(ngk)) ALLOCATE ( ngk(nks))
+     
+     ALLOCATE ( gk(npwx) )
+     igk_k(:,:) = 0
+     !
+     ! ... The following loop must NOT be called more than once in a run
+     ! ... or else there will be problems with variable-cell calculations
+     !
+     DO ik = 1, nks
+        CALL gk_sort( xk(1,ik), ngm, g, gcutw, ngk(ik), igk_k(1,ik), gk )
+     END DO
+     DEALLOCATE ( gk )
+     !
+   END SUBROUTINE init_igk
+   !
+   SUBROUTINE deallocate_igk ( ) 
+   IF ( ALLOCATED( ngk ) )        DEALLOCATE( ngk )
+   IF ( ALLOCATED( igk_k ) )      DEALLOCATE( igk_k )
+   END SUBROUTINE deallocate_igk
+ 
  END MODULE klist
  !
  MODULE lsda_mod
    !
    ! ... The variables needed for the lsda calculation
*** /home/giannozz/espresso-old/espresso-5.4.0/PW/src/read_file.f90	2016-04-24 23:19:18.000000000 +0200
--- PW/src/read_file.f90	2016-04-29 16:33:47.725714792 +0200
***************
*** 16,22 ****
    USE buffers,              ONLY : open_buffer, close_buffer
    USE wvfct,                ONLY : nbnd, npwx
    USE noncollin_module,     ONLY : npol
-   USE klist,                ONLY : nks
    USE paw_variables,        ONLY : okpaw, ddd_PAW
    USE paw_onecenter,        ONLY : paw_potential
    USE uspp,                 ONLY : becsum
--- 16,21 ----
***************
*** 27,32 ****
--- 26,34 ----
    USE ldaU,                 ONLY : lda_plus_u, U_projection
    USE pw_restart,           ONLY : pw_readfile
    USE control_flags,        ONLY : io_level
+   USE klist,                ONLY : init_igk
+   USE gvect,                ONLY : ngm, g
+   USE gvecw,                ONLY : gcutw
    !
    IMPLICIT NONE 
    INTEGER :: ierr
***************
*** 48,53 ****
--- 50,60 ----
    io_level = 1
    CALL open_buffer ( iunwfc, 'wfc', nwordwfc, io_level, exst )
    !
+   ! ... Allocate and compute k+G indices and number of plane waves
+   ! ... FIXME: should be read from file, not re-computed
+   !
+   CALL init_igk ( npwx, ngm, g, gcutw ) 
+   !
    ! ... Read orbitals, write them in 'distributed' form to iunwfc
    !
    CALL pw_readfile( 'wave', ierr )