!
! Copyright (C) 2005-2014 Quantum ESPRESSO group
! 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 .
!
!----------------------------------------------------------------------------
! TB
! included monopole related variables, search for 'TB'
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
Module pw_restart
    !----------------------------------------------------------------------------
    !
    ! ... this module contains methods to read and write data produced by PWscf
    !
    ! ... originally written by Carlo Sbraccia  (2005)
    !
#if ! defined(__XSD)
!
    Use iotk_module
    !
    Use qexml_module, Only: qexml_init, qexml_openfile, qexml_closefile, qexml_write_header, qexml_write_control, qexml_write_cell, qexml_write_moving_cell, qexml_write_ions, qexml_write_symmetry, qexml_write_efield, qexml_write_planewaves, qexml_write_spin, qexml_write_magnetization, qexml_write_xc, qexml_write_exx, qexml_write_occ, qexml_write_bz, qexml_write_para, qexml_write_bands_info, qexml_write_bands_pw, qexml_write_esm, qexml_wfc_filename, default_fmt_version => qexml_default_version, qexml_kpoint_dirname, qexml_read_header, qexml_read_cell, qexml_read_moving_cell, qexml_read_planewaves, qexml_read_ions, qexml_read_spin, qexml_read_magnetization, qexml_read_xc, qexml_read_occ, qexml_read_bz, qexml_read_bands_info, qexml_read_bands_pw, qexml_read_symmetry, qexml_read_efield, qexml_read_para, qexml_read_exx, qexml_read_esm
    !
    Use xml_io_base, Only: rho_binary, read_wfc, write_wfc, create_directory
    !
    !
    Use kinds, Only: DP
    Use constants, Only: e2, PI
    !
    Use io_files, Only: tmp_dir, prefix, iunpun, xmlpun, delete_if_present, qexml_version, qexml_version_init, pseudo_dir
    !
    Use io_global, Only: ionode, ionode_id
    Use mp_images, Only: intra_image_comm
    Use mp_pools, Only: my_pool_id
    Use mp_bands, Only: intra_bgrp_comm
    Use mp, Only: mp_bcast, mp_sum, mp_max
    Use parser, Only: version_compare
    !
    !
    Implicit None
    !
    Character (Len=256), External :: trimcheck
    !
    Save
    !
    Private
    !
    Public :: pw_readfile, pw_writefile
    Public :: gk_l2gmap, gk_l2gmap_kdip
    !
    Integer, Private :: iunout
    !
    Logical :: lcell_read = .False., lpw_read = .False., lions_read = .False., lspin_read = .False., lstarting_mag_read = .False., lxc_read = .False., locc_read = .False., lbz_read = .False., lbs_read = .False., lefield_read = .False., lwfc_read = .False., lsymm_read = .False.
    !
    !
Contains
!
!------------------------------------------------------------------------
    Subroutine pw_writefile (what)
    !------------------------------------------------------------------------
    !
        Use control_flags, Only: twfcollect, conv_ions, lscf, lkpoint_dir, gamma_only, tqr, tq_smoothing, tbeta_smoothing, noinv, do_makov_payne, smallmem, llondon, lxdm, ts_vdw
        Use realus, Only: real_space
        Use global_version, Only: version_number
        Use cell_base, Only: at, bg, alat, tpiba, tpiba2, ibrav, celldm
        Use gvect, Only: ig_l2g
        Use ions_base, Only: nsp, ityp, atm, nat, tau, if_pos
        Use noncollin_module, Only: noncolin, npol
        Use io_files, Only: nwordwfc, iunwfc, psfile
        Use buffers, Only: get_buffer
        Use wavefunctions_module, Only: evc
        Use klist, Only: nks, nkstot, xk, ngk, igk_k, wk, qnorm, lgauss, ngauss, degauss, nelec, two_fermi_energies, nelup, neldw
        Use start_k, Only: nk1, nk2, nk3, k1, k2, k3, nks_start, xk_start, wk_start
        Use ktetra, Only: ntetra, tetra
        Use klist, Only: ltetra
        Use gvect, Only: ngm, ngm_g, g, mill
        Use fft_base, Only: dfftp
        Use basis, Only: natomwfc
        Use gvecs, Only: ngms_g, dual
        Use fft_base, Only: dffts
        Use wvfct, Only: npw, npwx, et, wg, nbnd
        Use gvecw, Only: ecutwfc
        Use ener, Only: ef, ef_up, ef_dw
        Use fixed_occ, Only: tfixed_occ, f_inp
        Use ldaU, Only: lda_plus_u, lda_plus_u_kind, U_projection, Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_J, Hubbard_alpha, Hubbard_J0, Hubbard_beta
        Use spin_orb, Only: lspinorb, domag, lforcet
        Use symm_base, Only: nrot, nsym, invsym, s, ft, irt, t_rev, sname, time_reversal, no_t_rev
        Use lsda_mod, Only: nspin, isk, lsda, starting_magnetization
        Use noncollin_module, Only: angle1, angle2, i_cons, mcons, bfield, lambda
        Use ions_base, Only: amass
        Use funct, Only: get_dft_name, get_inlc
        Use kernel_table, Only: vdw_table_name
        Use scf, Only: rho
        Use extfield, Only: tefield, dipfield, edir, emaxpos, eopreg, eamp, monopole, zmon, block, block_1, block_2, block_height, relaxz
        Use io_rho_xml, Only: write_rho
        Use mp_world, Only: nproc
        Use mp_images, Only: nproc_image
        Use mp_pools, Only: kunit, nproc_pool, me_pool, root_pool, intra_pool_comm, inter_pool_comm
        Use mp_bands, Only: nproc_bgrp, me_bgrp, root_bgrp, intra_bgrp_comm, nbgrp, ntask_groups
        Use mp_diag, Only: nproc_ortho
        Use funct, Only: get_exx_fraction, dft_is_hybrid, get_gau_parameter, get_screening_parameter, exx_is_active
        Use exx, Only: x_gamma_extrapolation, nq1, nq2, nq3, exxdiv_treatment, yukawa, ecutvcut, ecutfock
        Use cellmd, Only: lmovecell, cell_factor
        Use martyna_tuckerman, Only: do_comp_mt
        Use esm, Only: do_comp_esm, esm_nfit, esm_efield, esm_w, esm_a, esm_bc
        Use acfdt_ener, Only: acfdt_in_pw
        Use london_module, Only: scal6, lon_rcut
        Use tsvdw_module, Only: vdw_isolated
        !
        !
        Implicit None
        !
        Character (Len=*), Intent (In) :: what
        !
        Character (Len=20) :: dft_name
        Character (Len=256) :: dirname, filename
        Integer :: i, ig, ik, ngg, ierr, ipol, num_k_points
        Integer :: nkl, nkr, npwx_g
        Integer :: ike, iks, npw_g, ispin, inlc
        Integer, External :: global_kpoint_index
        Integer, Allocatable :: ngk_g (:)
        Integer, Allocatable :: igk_l2g (:, :), igk_l2g_kdip (:, :), mill_g (:, :)
        Logical :: lwfc, lrho, lxsd
        Character (iotk_attlenx) :: attr
        !
        !
        Select Case (what)
        Case ("all")
        !
        ! ... do not overwrite the scf charge density with a non-scf one
        ! ... (except in the 'force theorem' calculation of MAE where the
        ! ...  charge density differs from the one read from disk)
        !
            lrho = lscf .Or. lforcet
            lwfc = twfcollect
            !
        Case ("config")
        !
        ! ... write just the xml data file, not the charge density and the wavefunctions
        !
            lwfc = .False.
            lrho = .False.
            !
        Case Default
        !
            Call errore ('pw_writefile', 'unexpected case: '// TRIM(what), 1)
            !
        End Select
        !
        If (ionode) Then
        !
        ! ... look for an empty unit (only ionode needs it)
        !
            Call iotk_free_unit (iunout, ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        Call errore ('pw_writefile ', 'no free units to write wavefunctions', ierr)
        !
        dirname = TRIM (tmp_dir) // TRIM (prefix) // '.save'
        !
        ! ... create the main restart directory
        !
        Call create_directory (dirname)
        !
        ! ... create the k-points subdirectories
        !
        If (nspin == 2) Then
            num_k_points = nkstot / 2
        Else
            num_k_points = nkstot
        End If
        !
        If (lkpoint_dir) Then
        !
            Do i = 1, num_k_points
            !
                Call create_directory (qexml_kpoint_dirname(dirname, i))
                !
            End Do
            !
        End If
        !
        iks = global_kpoint_index (nkstot, 1)
        ike = iks + nks - 1
        !
        ! ... find out the global number of G vectors: ngm_g
        !
        ngm_g = ngm
        !
        Call mp_sum (ngm_g, intra_bgrp_comm)
        !
        ! ... collect all G-vectors across processors within the pools
        !
        Allocate (mill_g(3, ngm_g))
        !
        mill_g = 0
        !
        Do ig = 1, ngm
        !
            mill_g (1, ig_l2g(ig)) = mill (1, ig)
            mill_g (2, ig_l2g(ig)) = mill (2, ig)
            mill_g (3, ig_l2g(ig)) = mill (3, ig)
            !
        End Do
        !
        Call mp_sum (mill_g, intra_bgrp_comm)
        !
        ! ... build the igk_l2g array, yielding the correspondence between
        ! ... the local k+G index and the global G index - see also ig_l2g
        ! ... igk_l2g is build from arrays igk, previously stored in hinit0
        ! ... Beware: for variable-cell case, one has to use starting G and
        ! ... k+G vectors
        !
        Allocate (igk_l2g(npwx, nks))
        igk_l2g = 0
        !
        Do ik = 1, nks
            npw = ngk (ik)
            Call gk_l2gmap (ngm, ig_l2g(1), npw, igk_k(1, ik), igk_l2g(1, ik))
        End Do
        !
        ! ... compute the global number of G+k vectors for each k point
        !
        Allocate (ngk_g(nkstot))
        !
        ngk_g = 0
        ngk_g (iks:ike) = ngk (1:nks)
        !
        Call mp_sum (ngk_g, inter_pool_comm)
        Call mp_sum (ngk_g, intra_pool_comm)
        !
        ngk_g = ngk_g / nbgrp
        !
        ! ... compute the maximum G vector index among all G+k an processors
        !
        npw_g = MAXVAL (igk_l2g(:, :))
        !
        Call mp_max (npw_g, inter_pool_comm)
        Call mp_max (npw_g, intra_pool_comm)
        !
        ! ... compute the maximum number of G vector among all k points
        !
        npwx_g = MAXVAL (ngk_g(1:nkstot))
        !
        ! ... define a further l2g map to write gkvectors and wfc coherently
        !
        Allocate (igk_l2g_kdip(npwx_g, nks))
        !
        igk_l2g_kdip = 0
        !
        Do ik = iks, ike
        !
            Call gk_l2gmap_kdip (npw_g, ngk_g(ik), ngk(ik-iks+1), igk_l2g(1, ik-iks+1), igk_l2g_kdip(1, ik-iks+1))
        End Do
        !
        If (ionode) Then
        !
        ! ... open XML descriptor
        !
            Call qexml_init (iunpun)
            Call qexml_openfile (TRIM(dirname)//'/'//TRIM(xmlpun), 'write', BINARY=.False., ierr=ierr)
            !
            If ( .Not. (lkpoint_dir)) Call iotk_open_write (iunout, FILE=TRIM(dirname)//'/'//TRIM(xmlpun)//'.eig', BINARY=.False., ierr=ierr)
            !
        End If
        !
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        Call errore ('pw_writefile ', 'cannot open restart file for writing', ierr)
        !
        If (ionode) Then
        !
        ! ... here we start writing the punch-file
        !
        !-------------------------------------------------------------------------------
        ! ... HEADER
        !-------------------------------------------------------------------------------
        !
            Call qexml_write_header ("PWSCF", TRIM(version_number))
            !
            !-------------------------------------------------------------------------------
            ! ... CONTROL
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_control (PP_CHECK_FLAG=conv_ions, lkpoint_dir=lkpoint_dir, Q_REAL_SPACE=tqr, tq_smoothing=tq_smoothing, BETA_REAL_SPACE=real_space, tbeta_smoothing=tbeta_smoothing)
            !
            !-------------------------------------------------------------------------------
            ! ... CELL
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_cell (ibrav, celldm, alat, at(:, 1), at(:, 2), at(:, 3), bg(:, 1), bg(:, 2), bg(:, 3), "Bohr", "Bohr",  "2 pi / a", do_makov_payne, do_comp_mt, do_comp_esm)
            !
            If (lmovecell) Call qexml_write_moving_cell (lmovecell, cell_factor)
            !
            !-------------------------------------------------------------------------------
            ! ... IONS
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_ions (nsp, nat, atm, ityp, psfile, pseudo_dir, amass, 'a.m.u.', tau, 'Bohr', if_pos, dirname, alat)
            !
            !-------------------------------------------------------------------------------
            ! ... SYMMETRIES
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_symmetry (ibrav, nrot, nsym, invsym, noinv, time_reversal, no_t_rev, ft, s, sname, "Crystal", irt, nat, t_rev)
            !
            !-------------------------------------------------------------------------------
            ! ... ELECTRIC FIELD
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_efield (tefield, dipfield, edir, emaxpos, eopreg, eamp, monopole, zmon, relaxz, block, block_1, block_2, block_height)
            !
            !
            !-------------------------------------------------------------------------------
            ! ... PLANE_WAVES
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_planewaves (ecutwfc/e2, ecutwfc*dual/e2, npwx_g, gamma_only, dfftp%nr1, dfftp%nr2, dfftp%nr3, ngm_g, dffts%nr1, dffts%nr2, dffts%nr3, ngms_g, dfftp%nr1, dfftp%nr2, dfftp%nr3, mill_g, lwfc, 'Hartree')
            !
            !-------------------------------------------------------------------------------
            ! ... SPIN
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_spin (lsda, noncolin, npol, lspinorb, domag)
            !
            Call qexml_write_magnetization (starting_magnetization, angle1*180.0_DP/PI, angle2*180.0_DP/PI, nsp, two_fermi_energies, i_cons, mcons, bfield, ef_up/e2, ef_dw/e2, nelup, neldw, lambda, 'Hartree')
            !
            !-------------------------------------------------------------------------------
            ! ... EXCHANGE_CORRELATION
            !-------------------------------------------------------------------------------
            !
            dft_name = get_dft_name ()
            inlc = get_inlc ()
            !
            Call qexml_write_xc (DFT=dft_name, nsp=nsp, lda_plus_u=lda_plus_u, lda_plus_u_kind=lda_plus_u_kind, U_projection=U_projection, Hubbard_lmax=Hubbard_lmax, Hubbard_l=Hubbard_l, Hubbard_U=Hubbard_U, Hubbard_J=Hubbard_J, Hubbard_J0=Hubbard_J0, Hubbard_beta=Hubbard_beta, Hubbard_alpha=Hubbard_alpha, inlc=inlc, vdw_table_name=vdw_table_name, pseudo_dir=pseudo_dir, dirname=dirname, acfdt_in_pw=acfdt_in_pw, llondon=llondon, LONDON_S6=scal6, LONDON_RCUT=lon_rcut, lxdm=lxdm, ts_vdw=ts_vdw, vdw_isolated=vdw_isolated)
            !
            !
            If (dft_is_hybrid()) Call qexml_write_exx (x_gamma_extrapolation, nq1, nq2, nq3, exxdiv_treatment, yukawa, ecutvcut, get_exx_fraction(), get_gau_parameter(), get_screening_parameter(), exx_is_active(), ecutfock)
            !
            !-------------------------------------------------------------------------------
            ! ... ESM
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_esm (esm_nfit, esm_efield, esm_w, esm_a, esm_bc)
            !
            !-------------------------------------------------------------------------------
            ! ... OCCUPATIONS
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_occ (lgauss=lgauss, ngauss=ngauss, degauss=degauss/e2, DEGAUSS_UNITS='Hartree', ltetra=ltetra, ntetra=ntetra, tetra=tetra, tfixed_occ=tfixed_occ, lsda=lsda, NSTATES_UP=nbnd, NSTATES_DW=nbnd, INPUT_OCC=f_inp)
            !
            !-------------------------------------------------------------------------------
            ! ... BRILLOUIN_ZONE
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_bz (num_k_points, xk, wk, k1, k2, k3, nk1, nk2, nk3, '2 pi / a', qnorm, nks_start, xk_start, wk_start)
            !
            !-------------------------------------------------------------------------------
            ! ... PARALLELISM
            !-------------------------------------------------------------------------------
            !
            !
            Call qexml_write_para (kunit, nproc, nproc_pool, nproc_image, ntask_groups, nproc_bgrp, nproc_ortho)
            !
            !-------------------------------------------------------------------------------
            ! ... CHARGE DENSITY
            !-------------------------------------------------------------------------------
            !
            !
            filename = "./charge-density.dat"
            If ( .Not. rho_binary) filename = "./charge-density.xml"
            !
            Call iotk_link (iunpun, "CHARGE-DENSITY", TRIM(filename), CREATE=.False., BINARY=.True.)
            !
            !-------------------------------------------------------------------------------
            ! ... BAND_STRUCTURE_INFO
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_bands_info (num_k_points, natomwfc, nbnd, nbnd, nbnd, nspin, nelec, Nint(nelup), Nint(neldw), "Hartree", "2 pi / a", ef=ef/e2, two_fermi_energies=two_fermi_energies, ef_up=ef_up/e2, ef_down=ef_dw/e2, noncolin=noncolin)
            !
            !-------------------------------------------------------------------------------
            ! ... EIGENVALUES
            !-------------------------------------------------------------------------------
            !
            Call qexml_write_bands_pw (nbnd, num_k_points, nspin, xk, wk, wg, et/e2, "Hartree", lkpoint_dir, iunout, dirname)
            !
            !
            If ( .Not. lkpoint_dir) Call iotk_close_write (iunout)
            !
            !-------------------------------------------------------------------------------
            ! ... EIGENVECTORS
            !-------------------------------------------------------------------------------
            !
            Call iotk_write_begin (iunpun, "EIGENVECTORS")
            !
            Call iotk_write_dat (iunpun, "MAX_NUMBER_OF_GK-VECTORS", npwx_g)
            !
        End If
        !
        k_points_loop2: Do ik = 1, num_k_points
        !
            If (ionode) Then
            !
                Call iotk_write_begin (iunpun, "K-POINT"// TRIM(iotk_index(ik)))
                !
                ! ... G+K vectors
                !
                Call iotk_write_dat (iunpun, "NUMBER_OF_GK-VECTORS", ngk_g(ik))
                !
                If (lwfc) Then
                !
                    filename = qexml_wfc_filename (".", 'gkvectors', ik, DIR=lkpoint_dir)
                    !
                    Call iotk_link (iunpun, "GK-VECTORS", filename, CREATE=.False., BINARY=.True.)
                    !
                    filename = qexml_wfc_filename (dirname, 'gkvectors', ik, DIR=lkpoint_dir)
                End If
                !
            End If
            !
            If (lwfc) Then
            !
                If ( .Not. smallmem) Call write_gk (iunout, ik, filename)
                !
                Call write_this_wfc (iunout, ik)
                !
            End If
            !
            If (ionode) Then
            !
                Call iotk_write_end (iunpun, "K-POINT"// TRIM(iotk_index(ik)))
                !
            End If
            !
        End Do k_points_loop2
        !
        If (ionode) Then
        !
            Call iotk_write_end (iunpun, "EIGENVECTORS")
            !
            Call qexml_closefile ('write', ierr=ierr)
            !
            !
            Call delete_if_present (TRIM(dirname)//'/'//TRIM(xmlpun)//'.bck')
            !
        End If
        !
        Deallocate (igk_l2g)
        Deallocate (igk_l2g_kdip)
        !
        !-------------------------------------------------------------------------------
        ! ... CHARGE-DENSITY FILES
        !-------------------------------------------------------------------------------
        !
        ! ... also writes rho%ns if lda+U and rho%bec if PAW
        !
        If (lrho) Call write_rho (rho, nspin)
        !-------------------------------------------------------------------------------
        ! ... END RESTART SECTIONS
        !-------------------------------------------------------------------------------
        !
        Deallocate (mill_g)
        Deallocate (ngk_g)
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        Call errore ('pw_writefile ', 'cannot save history', ierr)
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        Call errore ('pw_writefile ', 'cannot save history', ierr)
        !
        Return
        !
    Contains
    !
    !--------------------------------------------------------------------
        Subroutine write_gk (iun, ik, filename)
        !--------------------------------------------------------------------
        !
#if defined __HDF5
            Use hdf5_qe, Only: prepare_for_writing_final, write_gkhdf5, gk_hdf5_write, h5fclose_f
            Use io_files, Only: tmp_dir
#endif
!
            Implicit None
            !
            Integer, Intent (In) :: iun, ik
            Character (Len=256), Intent (In) :: filename
            !
            Integer, Allocatable :: igwk (:, :)
            Integer, Allocatable :: itmp (:)
            Character (Len=256) :: filename_hdf5
            Integer :: ierr
            !
            !
            Allocate (igwk(npwx_g, nkstot))
            !
            igwk (:, ik) = 0
            !
            Allocate (itmp(npw_g))
            !
            itmp = 0
            !
            If (ik >= iks .And. ik <= ike) Then
            !
                Do ig = 1, ngk (ik-iks+1)
                !
                    itmp (igk_l2g(ig, ik-iks+1)) = igk_l2g (ig, ik-iks+1)
                    !
                End Do
                !
            End If
            !
            Call mp_sum (itmp, inter_pool_comm)
            Call mp_sum (itmp, intra_pool_comm)
            !
            ngg = 0
            !
            Do ig = 1, npw_g
            !
                If (itmp(ig) == ig) Then
                !
                    ngg = ngg + 1
                    !
                    igwk (ngg, ik) = ig
                    !
                End If
                !
            End Do
            !
            Deallocate (itmp)
            !
            If (ionode) Then
            !
#if defined __HDF5
                filename_hdf5 = TRIM (tmp_dir) // "gk.hdf5"
                Call prepare_for_writing_final (gk_hdf5_write, inter_pool_comm, filename_hdf5, ik)
                Call write_gkhdf5 (gk_hdf5_write, xk(:, ik), igwk(1:ngk_g(ik), ik), mill_g(1:3, igwk(1:ngk_g(ik), ik)), ik)
                Call h5fclose_f (gk_hdf5_write%file_id, ierr)
#else
!
                Call iotk_open_write (iun, FILE=TRIM(filename), ROOT="GK-VECTORS", BINARY=.True.)
                !
                Call iotk_write_dat (iun, "NUMBER_OF_GK-VECTORS", ngk_g(ik))
                Call iotk_write_dat (iun, "MAX_NUMBER_OF_GK-VECTORS", npwx_g)
                Call iotk_write_dat (iun, "GAMMA_ONLY", gamma_only)
                !
                Call iotk_write_attr (attr, "UNITS", "2 pi / a", FIRST=.True.)
                Call iotk_write_dat (iun, "K-POINT_COORDS", xk(:, ik), attr=attr)
                !
                Call iotk_write_dat (iun, "INDEX", igwk(1:ngk_g(ik), ik))
                Call iotk_write_dat (iun, "GRID", mill_g(1:3, igwk(1:ngk_g(ik), ik)), COLUMNS=3)
                !
                Call iotk_close_write (iun)
                !
#endif
            End If
            !
            Deallocate (igwk)
            !
        End Subroutine write_gk
        !
        !
        !--------------------------------------------------------------------
        Subroutine write_this_wfc (iun, ik)
        !--------------------------------------------------------------------
        !
            Implicit None
            !
            Integer, Intent (In) :: iun, ik
            Character (Len=256) :: filename
            Integer :: ispin, ik_eff
            !
            ! ... wavefunctions - do not read if already in memory (nsk==1)
            ! ...                 read only if on this pool (iks <= ik <= ike )
            !
            If ((nks > 1) .And. (ik >= iks) .And. (ik <= ike)) Then
                Call get_buffer (evc, nwordwfc, iunwfc, (ik-iks+1))
            End If
            !
            If (nspin == 2) Then
            !
            ! ... beware: with pools, isk(ik) has the correct value for
            ! ... all k-points only on first pool (ionode proc is ok)
            !
                ispin = isk (ik)
                !
                If (ionode) Then
                !
                    filename = qexml_wfc_filename (".", 'evc', ik, ispin, DIR=lkpoint_dir)
                    !
                    Call iotk_link (iunpun, "WFC"// TRIM(iotk_index(ispin)), filename, CREATE=.False., BINARY=.True.)
                    !
                    filename = qexml_wfc_filename (dirname, 'evc', ik, ispin, DIR=lkpoint_dir)
                    !
                End If
                !
                Call write_wfc (iunout, ik, nkstot, kunit, ispin, nspin, evc, npw_g, gamma_only, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, 1.D0, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                !
                ik_eff = ik + num_k_points
                !
                ispin = isk (ik_eff)
                !
                ! ... LSDA: now read minority wavefunctions (if not already
                ! ... in memory and if they are on this pool)
                !
                If ((nks > 1) .And. (ik_eff >= iks) .And. (ik_eff <= ike)) Then
                !
                    Call get_buffer (evc, nwordwfc, iunwfc, (ik_eff-iks+1))
                    !
                End If
                !
                If (ionode) Then
                !
                    filename = qexml_wfc_filename (".", 'evc', ik, ispin, DIR=lkpoint_dir)
                    !
                    Call iotk_link (iunpun, "WFC"// TRIM(iotk_index(ispin)), filename, CREATE=.False., BINARY=.True.)
                    !
                    filename = qexml_wfc_filename (dirname, 'evc', ik, ispin, DIR=lkpoint_dir)
                    !
                End If
                !
                Call write_wfc (iunout, ik_eff, nkstot, kunit, ispin, nspin, evc, npw_g, gamma_only, nbnd, igk_l2g_kdip(:, ik_eff-iks+1), ngk(ik_eff-iks+1), filename, 1.D0, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                !
            Else
            !
                If (noncolin) Then
                !
                    Do ipol = 1, npol
                    !
                        If (ionode) Then
                        !
                            filename = qexml_wfc_filename (".", 'evc', ik, ipol, DIR=lkpoint_dir)
                            !
                            Call iotk_link (iunpun, "WFC"// TRIM(iotk_index(ipol)), filename, CREATE= .False., BINARY=.True.)
                            !
                            filename = qexml_wfc_filename (dirname, 'evc', ik, ipol, DIR=lkpoint_dir)
                            !
                        End If
                        !
                        ! TEMP  spin-up and spin-down spinor components are written
                        ! TEMP  to different files, like in LSDA - not a smart way
                        !
                        nkl = (ipol-1) * npwx + 1
                        nkr = ipol * npwx
                        Call write_wfc (iunout, ik, nkstot, kunit, ipol, npol, evc(nkl:nkr, :), npw_g, gamma_only, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, 1.D0, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                        !
                    End Do
                    !
                Else
                !
                    ispin = 1
                    !
                    If (ionode) Then
                    !
                        filename = qexml_wfc_filename (".", 'evc', ik, DIR=lkpoint_dir)
                        !
                        Call iotk_link (iunpun, "WFC", filename, CREATE=.False., BINARY=.True.)
                        !
                        filename = qexml_wfc_filename (dirname, 'evc', ik, DIR=lkpoint_dir)
                        !
                    End If
                    !
                    Call write_wfc (iunout, ik, nkstot, kunit, ispin, nspin, evc, npw_g, gamma_only, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, 1.D0, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                    !
                End If
                !
            End If
            !
        End Subroutine write_this_wfc
        !
    End Subroutine pw_writefile
    !
    !------------------------------------------------------------------------
    Subroutine pw_readfile (what, ierr)
    !------------------------------------------------------------------------
    !
        Use io_rho_xml, Only: read_rho
        Use scf, Only: rho
        Use lsda_mod, Only: nspin
        Use mp_bands, Only: intra_bgrp_comm
        Use mp, Only: mp_sum
        !
        Implicit None
        !
        Character (Len=*), Intent (In) :: what
        Integer, Intent (Out) :: ierr
        !
        Character (Len=256) :: dirname
        Character (Len=80) :: errmsg
        Logical :: lcell, lpw, lions, lspin, linit_mag, lxc, locc, lbz, lbs, lwfc, lheader, lsymm, lrho, lefield, ldim, lef, lexx, lesm
        !
        Integer :: tmp
        !
        ierr = 0
        !
        dirname = TRIM (tmp_dir) // TRIM (prefix) // '.save'
        !
        ! ... look for an empty unit
        !
        Call iotk_free_unit (iunout, ierr)
        !
        Call errore ('pw_readfile', 'no free units to read wavefunctions', ierr)
        !
        lheader = .Not. qexml_version_init
        !
        ldim = .False.
        lcell = .False.
        lpw = .False.
        lions = .False.
        lspin = .False.
        linit_mag = .False.
        lxc = .False.
        locc = .False.
        lbz = .False.
        lbs = .False.
        lwfc = .False.
        lsymm = .False.
        lrho = .False.
        lefield = .False.
        lef = .False.
        lexx = .False.
        lesm = .False.
        !
        Select Case (what)
        Case ('header')
        !
            lheader = .True.
            !
        Case ('dim')
        !
            ldim = .True.
            lbz = .True.
            !
        Case ('pseudo')
        !
            lions = .True.
            !
        Case ('config')
        !
            lcell = .True.
            lions = .True.
            !
        Case ('wave')
        !
            lpw = .True.
            lwfc = .True.
            !
        Case ('nowavenobs')
        !
            lcell = .True.
            lpw = .True.
            lions = .True.
            lspin = .True.
            linit_mag = .True.
            lxc = .True.
            locc = .True.
            lbz = .True.
            lsymm = .True.
            lefield = .True.
            !
        Case ('nowave')
        !
            lcell = .True.
            lpw = .True.
            lions = .True.
            lspin = .True.
            linit_mag = .True.
            lxc = .True.
            locc = .True.
            lbz = .True.
            lbs = .True.
            lsymm = .True.
            lefield = .True.
            !
        Case ('all')
        !
            lcell = .True.
            lpw = .True.
            lions = .True.
            lspin = .True.
            linit_mag = .True.
            lxc = .True.
            locc = .True.
            lbz = .True.
            lbs = .True.
            lwfc = .True.
            lsymm = .True.
            lefield = .True.
            lrho = .True.
            !
        Case ('reset')
        !
            lcell_read = .False.
            lpw_read = .False.
            lions_read = .False.
            lspin_read = .False.
            lstarting_mag_read = .False.
            lxc_read = .False.
            locc_read = .False.
            lbz_read = .False.
            lbs_read = .False.
            lwfc_read = .False.
            lsymm_read = .False.
            lefield_read = .False.
            !
        Case ('ef')
        !
            lef = .True.
            !
        Case ('exx')
        !
            lexx = .True.
            !
        Case ('esm')
        !
            lesm = .True.
            !
        Case Default
        !
            Call errore ('pw_readfile', 'unknown case '//TRIM(what), 1)
            !
        End Select
        !
        If ( .Not. lheader .And. .Not. qexml_version_init) Call errore ('pw_readfile', 'qexml version not set', 71)
        !
        If (ionode) Then
        !
            Call qexml_init (iunpun)
            Call qexml_openfile (TRIM(dirname)//'/'//TRIM(xmlpun), 'read', BINARY=.False., ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr /= 0) Then
            errmsg = 'error opening xml data file'
            Go To 100
        End If
        !
        If (lheader) Then
        !
            Call read_header (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading header of xml data file'
                Go To 100
            End If
            !
        End If
        !
        If (ldim) Then
        !
            Call read_dim (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading dimensions in xml data file'
                Go To 100
            End If
            !
        End If
        !
        If (lcell) Then
        !
            Call read_cell (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading cell info in xml data file'
                Go To 100
            End If
            !
        End If
        If (lpw) Then
        !
            Call read_planewaves (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading plane-wave info in xml data file'
                Go To 100
            End If
            !
        End If
        If (lions) Then
        !
            Call read_ions (dirname, ierr)
            If (ierr > 0) Then
                errmsg = 'error reading info on ions in xml data file'
                Go To 100
            End If
            !
        End If
        If (lspin) Then
        !
            Call read_spin (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading spin in xml data file'
                Go To 100
            End If
            !
        End If
        If (linit_mag) Then
        !
            Call read_magnetization (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading magnetization in xml data file'
                Go To 100
            End If
            !
        End If
        If (lxc) Then
        !
            Call read_xc (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading XC functional in xml data file'
                Go To 100
            End If
            !
        End If
        If (locc) Then
        !
            Call read_occupations (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading occupation numbers in xml data file'
                Go To 100
            End If
            !
        End If
        If (lbz) Then
        !
            Call read_brillouin_zone (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading Brillouin Zone in xml data file'
                Go To 100
            End If
            !
        End If
        If (lbs) Then
        !
            Call read_band_structure (dirname, ierr)
            If (ierr > 0) Then
                errmsg = 'error reading band structure in xml data file'
                Go To 100
            End If
            !
        End If
        If (lwfc) Then
        !
            Call read_wavefunctions (dirname, ierr)
            If (ierr > 0) Then
                errmsg = 'error reading wavefunctions in xml data file'
                Go To 100
            End If
            !
        End If
        If (lsymm) Then
        !
            Call read_symmetry (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading symmetry in xml data file'
                Go To 100
            End If
            !
        End If
        If (lefield) Then
        !
            Call read_efield (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading electric fields in xml data file'
                Go To 100
            End If
            !
        End If
        !
        If (lrho) Then
        !
        ! ... to read the charge-density we use the routine from io_rho_xml
        ! ... it also reads ns for ldaU and becsum for PAW
        !
            Call read_rho (rho, nspin)
            !
        End If
        !
        If (lef) Then
        !
            Call read_ef (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading Fermi energy and number of electrons in xml data file'
                Go To 100
            End If
            !
        End If
        If (lexx) Then
        !
            Call read_exx (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading hybrid functional in xml data file'
                Go To 100
            End If
            !
        End If
        If (lesm) Then
        !
            Call read_esm (ierr)
            If (ierr > 0) Then
                errmsg = 'error reading ESM restart data in xml data file'
                Go To 100
            End If
            !
        End If
        !
        If (ionode) Then
        !
            Call qexml_closefile ('read', ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Then
            errmsg = 'error closing xml data file'
            Go To 100
        End If
        !
        !
        Return
        !
        ! uncomment to continue execution after an error occurs
        ! 100 IF (ionode) THEN
        !        CALL qexml_closefile( 'read', IERR=tmp)
        !     ENDIF
        !     RETURN
        ! comment to continue execution after an error occurs
        !
        !
100      Call errore ('pw_readfile', TRIM(errmsg), ierr)
!
    End Subroutine pw_readfile
    !
    !------------------------------------------------------------------------
    Subroutine read_header (ierr)
    !------------------------------------------------------------------------
    !
    ! ... this routine reads the format version of the current xml datafile
    !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        ierr = 0
        !
        If (qexml_version_init) Return
        !
        If (ionode) Then
        !
            Call qexml_read_header (FORMAT_VERSION=qexml_version, ierr=ierr)
            !
            qexml_version_init = .True.
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr /= 0) Return
        !
        Call mp_bcast (qexml_version, ionode_id, intra_image_comm)
        Call mp_bcast (qexml_version_init, ionode_id, intra_image_comm)
        !
        !
    End Subroutine read_header
    !
    !------------------------------------------------------------------------
    Subroutine read_dim (ierr)
    !------------------------------------------------------------------------
    !
    ! ... this routine collects array dimensions from various sections
    ! ... plus with some other variables needed for array allocation
    !
        Use ions_base, Only: nat, nsp
        Use symm_base, Only: nsym
        Use gvect, Only: ngm_g, ecutrho
        Use fft_base, Only: dfftp
        Use gvecs, Only: ngms_g, dual
        Use fft_base, Only: dffts
        Use lsda_mod, Only: lsda
        Use noncollin_module, Only: noncolin
        Use ktetra, Only: ntetra
        Use klist, Only: nkstot, nelec
        Use wvfct, Only: nbnd, npwx
        Use gvecw, Only: ecutwfc
        Use control_flags, Only: gamma_only
        Use mp_pools, Only: kunit
        Use mp_global, Only: nproc_file, nproc_pool_file, nproc_image_file, ntask_groups_file, nproc_bgrp_file, nproc_ortho_file
        !
        Implicit None
        !
        !CHARACTER(LEN=*), INTENT(IN)  :: dirname
        Integer, Intent (Out) :: ierr
        !
        Integer :: npwx_
        Logical :: found, found2
        Character (iotk_attlenx) :: attr
        !
        !
        ! ... first the entire CELL section is read
        ! ...
        ierr = 0
        !
        Call read_cell (ierr)
        If (ierr /= 0) Go To 100
        !
        If (ionode) Then
        !
            Call qexml_read_ions (nat=nat, nsp=nsp, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            Call qexml_read_symmetry (nsym=nsym, found=found, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            If ( .Not. found) Then
            !
                nsym = 1
                !
            End If
            !
            Call qexml_read_planewaves (ecutwfc=ecutwfc, ecutrho=ecutrho, npwx=npwx_, gamma_only=gamma_only, nr1=dfftp%nr1, nr2=dfftp%nr2, nr3=dfftp%nr3, ngm=ngm_g, NR1S=dffts%nr1, NR2S=dffts%nr2, NR3S=dffts%nr3, NGMS=ngms_g, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            ecutwfc = ecutwfc * e2
            ecutrho = ecutrho * e2
            !
            dual = ecutrho / ecutwfc
            !
            Call qexml_read_spin (lsda=lsda, noncolin=noncolin, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            Call qexml_read_occ (ntetra=ntetra, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            Call qexml_read_bz (num_k_points=nkstot, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            If (lsda) nkstot = nkstot * 2
            !
            Call qexml_read_bands_info (nbnd=nbnd, nelec=nelec, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            Call qexml_read_para (kunit=kunit, nproc=nproc_file, nproc_pool=nproc_pool_file, nproc_image=nproc_image_file, ntask_groups=ntask_groups_file, nproc_bgrp=nproc_bgrp_file, nproc_ortho=nproc_ortho_file, found=found, ierr=ierr)
            If (ierr /= 0) Go To 100
            !
            If ( .Not. found) Then
            !
                kunit = 1
                nproc_file = 1
                nproc_pool_file = 1
                nproc_image_file = 1
                ntask_groups_file = 1
                nproc_bgrp_file = 1
                nproc_ortho_file = 1
                !
            End If
            !
        End If
        !
100      Call mp_bcast (ierr, ionode_id, intra_image_comm)
!
        If (ierr > 0) Return
        !
        Call mp_bcast (nat, ionode_id, intra_image_comm)
        Call mp_bcast (nsp, ionode_id, intra_image_comm)
        Call mp_bcast (nsym, ionode_id, intra_image_comm)
        Call mp_bcast (ecutwfc, ionode_id, intra_image_comm)
        Call mp_bcast (ecutrho, ionode_id, intra_image_comm)
        Call mp_bcast (dual, ionode_id, intra_image_comm)
        Call mp_bcast (npwx_, ionode_id, intra_image_comm)
        Call mp_bcast (gamma_only, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr1, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr2, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr3, ionode_id, intra_image_comm)
        Call mp_bcast (ngm_g, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr1, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr2, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr3, ionode_id, intra_image_comm)
        Call mp_bcast (ngms_g, ionode_id, intra_image_comm)
        Call mp_bcast (lsda, ionode_id, intra_image_comm)
        Call mp_bcast (noncolin, ionode_id, intra_image_comm)
        Call mp_bcast (ntetra, ionode_id, intra_image_comm)
        Call mp_bcast (nkstot, ionode_id, intra_image_comm)
        Call mp_bcast (nelec, ionode_id, intra_image_comm)
        Call mp_bcast (nbnd, ionode_id, intra_image_comm)
        Call mp_bcast (kunit, ionode_id, intra_image_comm)
        Call mp_bcast (nproc_file, ionode_id, intra_image_comm)
        Call mp_bcast (nproc_pool_file, ionode_id, intra_image_comm)
        Call mp_bcast (nproc_image_file, ionode_id, intra_image_comm)
        Call mp_bcast (ntask_groups_file, ionode_id, intra_image_comm)
        Call mp_bcast (nproc_bgrp_file, ionode_id, intra_image_comm)
        Call mp_bcast (nproc_ortho_file, ionode_id, intra_image_comm)
        !
        Return
        !
    End Subroutine read_dim
    !
    !--------------------------------------------------------------------------
    Subroutine read_cell (ierr)
    !------------------------------------------------------------------------
    !
        Use run_info, Only: title
        Use cell_base, Only: ibrav, alat, at, bg, celldm
        Use cell_base, Only: tpiba, tpiba2, omega
        Use cellmd, Only: lmovecell, cell_factor
        Use control_flags, Only: do_makov_payne
        Use martyna_tuckerman, Only: do_comp_mt
        Use esm, Only: do_comp_esm
        !
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        Character (Len=80) :: bravais_lattice, es_corr
        !
        !
        ierr = 0
        If (lcell_read) Return
        !
        If (ionode) Then
        !
            Call qexml_read_cell (bravais_lattice=bravais_lattice, celldm=celldm, alat=alat, A1=at(:, 1), A2=at(:, 2), A3=at(:, 3), es_corr=es_corr, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            Select Case (TRIM(es_corr))
            Case ("Makov-Payne")
                do_makov_payne = .True.
                do_comp_mt = .False.
                do_comp_esm = .False.
            Case ("Martyna-Tuckerman")
                do_makov_payne = .False.
                do_comp_mt = .True.
                do_comp_esm = .False.
            Case ("ESM")
                do_makov_payne = .False.
                do_comp_mt = .False.
                do_comp_esm = .True.
            Case ("None")
                do_makov_payne = .False.
                do_comp_mt = .False.
                do_comp_esm = .False.
            End Select
            !
            Select Case (TRIM(bravais_lattice))
            Case ("free")
                ibrav = 0
            Case ("cubic P (sc)")
                ibrav = 1
            Case ("cubic F (fcc)")
                ibrav = 2
            Case ("cubic I (bcc)")
                ibrav = 3
            Case ("Hexagonal and Trigonal P")
                ibrav = 4
            Case ("Trigonal R")
                ibrav = 5
            Case ("Tetragonal P (st)")
                ibrav = 6
            Case ("Tetragonal I (bct)")
                ibrav = 7
            Case ("Orthorhombic P")
                ibrav = 8
            Case ("Orthorhombic base-centered(bco)")
                ibrav = 9
            Case ("Orthorhombic face-centered")
                ibrav = 10
            Case ("Orthorhombic body-centered")
                ibrav = 11
            Case ("Monoclinic P")
                ibrav = 12
            Case ("Monoclinic base-centered")
                ibrav = 13
            Case ("Triclinic P")
                ibrav = 14
            Case Default
                ibrav = 0
            End Select
            !
            ! ... some internal variables
            !
            tpiba = 2.D0 * PI / alat
            tpiba2 = tpiba ** 2
            !
            ! ... to alat units
            !
            at (:, :) = at (:, :) / alat
            !
            Call volume (alat, at(1, 1), at(1, 2), at(1, 3), omega)
            !
            ! ... Generate the reciprocal lattice vectors
            !
            Call recips (at(1, 1), at(1, 2), at(1, 3), bg(1, 1), bg(1, 2), bg(1, 3))
            !
            Call qexml_read_moving_cell (lmovecell, cell_factor, ierr)
            !
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Return
        !
        Call mp_bcast (ibrav, ionode_id, intra_image_comm)
        Call mp_bcast (alat, ionode_id, intra_image_comm)
        Call mp_bcast (celldm, ionode_id, intra_image_comm)
        Call mp_bcast (tpiba, ionode_id, intra_image_comm)
        Call mp_bcast (tpiba2, ionode_id, intra_image_comm)
        Call mp_bcast (omega, ionode_id, intra_image_comm)
        Call mp_bcast (at, ionode_id, intra_image_comm)
        Call mp_bcast (bg, ionode_id, intra_image_comm)
        Call mp_bcast (do_makov_payne, ionode_id, intra_image_comm)
        Call mp_bcast (do_comp_mt, ionode_id, intra_image_comm)
        Call mp_bcast (do_comp_esm, ionode_id, intra_image_comm)
        Call mp_bcast (lmovecell, ionode_id, intra_image_comm)
        If (lmovecell) Then
            Call mp_bcast (cell_factor, ionode_id, intra_image_comm)
        Else
            cell_factor = 1.0_DP
        End If
        !
        title = ' '
        !
        lcell_read = .True.
        !
        Return
        !
    End Subroutine read_cell
    !
    !
    !------------------------------------------------------------------------
    Subroutine read_ions (dirname, ierr)
    !------------------------------------------------------------------------
    !
        Use ions_base, Only: nat, nsp, ityp, amass, atm, tau, if_pos
        Use cell_base, Only: alat
        Use io_files, Only: psfile, pseudo_dir, pseudo_dir_cur
        !
        Implicit None
        !
        Character (Len=*), Intent (In) :: dirname
        Integer, Intent (Out) :: ierr
        !
        Integer :: i
        Logical :: exst
        !
        ierr = 0
        If (lions_read) Return
        !
        If ( .Not. lcell_read) Call errore ('read_ions', 'read cell first', 1)
        !
        ! this is where PP files should be read from
        !
        pseudo_dir_cur = trimcheck (dirname)
        !
        If (ionode) Then
        !
            Call qexml_read_ions (nsp=nsp, nat=nat, atm=atm, ityp=ityp, psfile=psfile, amass=amass, tau=tau, if_pos=if_pos, pseudo_dir=pseudo_dir, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            Do i = 1, nat
            !
                tau (:, i) = tau (:, i) / alat
                !
            End Do
            !
        End If
        !
        Call mp_bcast (nat, ionode_id, intra_image_comm)
        Call mp_bcast (nsp, ionode_id, intra_image_comm)
        Call mp_bcast (atm, ionode_id, intra_image_comm)
        Call mp_bcast (amass, ionode_id, intra_image_comm)
        Call mp_bcast (psfile, ionode_id, intra_image_comm)
        Call mp_bcast (pseudo_dir, ionode_id, intra_image_comm)
        Call mp_bcast (ityp, ionode_id, intra_image_comm)
        Call mp_bcast (tau, ionode_id, intra_image_comm)
        Call mp_bcast (if_pos, ionode_id, intra_image_comm)
        !
        lions_read = .True.
        !
        Return
        !
    End Subroutine read_ions
    !
    !------------------------------------------------------------------------
    Subroutine read_symmetry (ierr)
    !------------------------------------------------------------------------
    !
        Use symm_base, Only: nrot, nsym, invsym, s, ft, ftau, irt, t_rev, sname, sr, invs, inverse_s, s_axis_to_cart, time_reversal, no_t_rev
        Use control_flags, Only: noinv
        Use fft_base, Only: dfftp
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        Character (iotk_attlenx) :: attr
        !
        Integer :: i
        Logical :: found
        !
        ierr = 0
        If (lsymm_read) Return
        !
        If ( .Not. lpw_read) Call errore ('read_symmetry', 'read planewaves first', 1)
        !
        If (ionode) Then
        !
            Call qexml_read_symmetry (nsym=nsym, nrot=nrot, invsym=invsym, noinv=noinv, time_reversal=time_reversal, no_t_rev=no_t_rev, TRASL=ft, s=s, sname=sname, t_rev=t_rev, irt=irt, found=found, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            If ( .Not. found) Then
            !
                nsym = 1
                s (:, :, nsym) = 0
                s (1, 1, nsym) = 1
                s (2, 2, nsym) = 1
                s (3, 3, nsym) = 1
                sr (:, :, nsym) = DBLE (s(:, :, nsym))
                ftau (:, nsym) = 0
                ft (:, nsym) = 0.0_DP
                sname (nsym) = 'identity'
                Do i = 1, SIZE (irt, 2)
                    irt (nsym, i) = i
                End Do
                invsym = .False.
                noinv = .False.
                t_rev (nsym) = 0
                invs (1) = 1
                time_reversal = .True.
                no_t_rev = .False.
                !
            Else
            !
                Do i = 1, nsym
                !
                    ftau (1, i) = Nint (ft(1, i)*DBLE(dfftp%nr1))
                    ftau (2, i) = Nint (ft(2, i)*DBLE(dfftp%nr2))
                    ftau (3, i) = Nint (ft(3, i)*DBLE(dfftp%nr3))
                    !
                End Do
                !
                ! indices of inverse operations and matrices in cartesian axis
                ! are not saved to disk (maybe they should), are recalculated here
                !
                Call inverse_s ()
                Call s_axis_to_cart ()
                !
            End If
            !
            !
        End If
        !
        Call mp_bcast (nsym, ionode_id, intra_image_comm)
        Call mp_bcast (nrot, ionode_id, intra_image_comm)
        Call mp_bcast (invsym, ionode_id, intra_image_comm)
        Call mp_bcast (noinv, ionode_id, intra_image_comm)
        Call mp_bcast (time_reversal, ionode_id, intra_image_comm)
        Call mp_bcast (no_t_rev, ionode_id, intra_image_comm)
        Call mp_bcast (s, ionode_id, intra_image_comm)
        Call mp_bcast (ftau, ionode_id, intra_image_comm)
        Call mp_bcast (ft, ionode_id, intra_image_comm)
        Call mp_bcast (sname, ionode_id, intra_image_comm)
        Call mp_bcast (irt, ionode_id, intra_image_comm)
        Call mp_bcast (t_rev, ionode_id, intra_image_comm)
        Call mp_bcast (invs, ionode_id, intra_image_comm)
        Call mp_bcast (sr, ionode_id, intra_image_comm)
        !
        lsymm_read = .True.
        !
        Return
        !
    End Subroutine read_symmetry
    !
    !------------------------------------------------------------------------
    Subroutine read_efield (ierr)
    !----------------------------------------------------------------------
    !
        Use extfield, Only: tefield, dipfield, edir, emaxpos, eopreg, eamp, monopole, zmon, relaxz, block, block_1, block_2, block_height
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        Logical :: found
        !
        ierr = 0
        If (lefield_read) Return
        !
        !
        If (ionode) Then
        !
            Call qexml_read_efield (tefield=tefield, dipfield=dipfield, edir=edir, emaxpos=emaxpos, eopreg=eopreg, eamp=eamp, monopole=monopole, zmon=zmon, relaxz=relaxz, block=block, block_1=block_1, block_2=block_2, block_height=block_height, found=found, ierr=ierr)
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If ((ionode) .And. ( .Not. found)) Then
        !
            tefield = .False.
            dipfield = .False.
            monopole = .False.
            !
        End If
        !
        Call mp_bcast (tefield, ionode_id, intra_image_comm)
        Call mp_bcast (dipfield, ionode_id, intra_image_comm)
        Call mp_bcast (edir, ionode_id, intra_image_comm)
        Call mp_bcast (emaxpos, ionode_id, intra_image_comm)
        Call mp_bcast (eopreg, ionode_id, intra_image_comm)
        Call mp_bcast (eamp, ionode_id, intra_image_comm)
        Call mp_bcast (monopole, ionode_id, intra_image_comm)
        Call mp_bcast (zmon, ionode_id, intra_image_comm)
        Call mp_bcast (relaxz, ionode_id, intra_image_comm)
        Call mp_bcast (block, ionode_id, intra_image_comm)
        Call mp_bcast (block_1, ionode_id, intra_image_comm)
        Call mp_bcast (block_2, ionode_id, intra_image_comm)
        Call mp_bcast (block_height, ionode_id, intra_image_comm)
        !
        lefield_read = .True.
        !
        Return
        !
    End Subroutine read_efield
    !
    !------------------------------------------------------------------------
    Subroutine read_planewaves (ierr)
    !------------------------------------------------------------------------
    !
        Use gvect, Only: ngm_g, ecutrho
        Use gvecs, Only: ngms_g, dual
        Use gvecw, Only: ecutwfc
        Use fft_base, Only: dfftp
        Use fft_base, Only: dffts
        Use wvfct, Only: npwx
        Use control_flags, Only: gamma_only
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        Integer :: npwx_
        !
        ierr = 0
        If (lpw_read) Return
        !
        !
        If (ionode) Call qexml_read_planewaves (ecutwfc=ecutwfc, ecutrho=ecutrho, npwx=npwx_, gamma_only=gamma_only, nr1=dfftp%nr1, nr2=dfftp%nr2, nr3=dfftp%nr3, ngm=ngm_g, NR1S=dffts%nr1, NR2S=dffts%nr2, NR3S=dffts%nr3, NGMS=ngms_g, ierr=ierr)
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            ecutwfc = ecutwfc * e2
            ecutrho = ecutrho * e2
            !
            dual = ecutrho / ecutwfc
            !
        End If
        !
        Call mp_bcast (ecutwfc, ionode_id, intra_image_comm)
        Call mp_bcast (ecutrho, ionode_id, intra_image_comm)
        Call mp_bcast (dual, ionode_id, intra_image_comm)
        Call mp_bcast (npwx_, ionode_id, intra_image_comm)
        Call mp_bcast (gamma_only, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr1, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr2, ionode_id, intra_image_comm)
        Call mp_bcast (dfftp%nr3, ionode_id, intra_image_comm)
        Call mp_bcast (ngm_g, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr1, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr2, ionode_id, intra_image_comm)
        Call mp_bcast (dffts%nr3, ionode_id, intra_image_comm)
        Call mp_bcast (ngms_g, ionode_id, intra_image_comm)
        !
        lpw_read = .True.
        !
        Return
        !
    End Subroutine read_planewaves
    !
    !------------------------------------------------------------------------
    Subroutine read_spin (ierr)
    !------------------------------------------------------------------------
    !
        Use spin_orb, Only: lspinorb, domag
        Use lsda_mod, Only: nspin, lsda
        Use noncollin_module, Only: noncolin, npol
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        Logical :: found
        !
        ierr = 0
        If (lspin_read) Return
        !
        If (ionode) Then
        !
            Call qexml_read_spin (lsda, noncolin, npol, lspinorb, domag, ierr)
            !
            If (lsda) Then
            !
                nspin = 2
                !
            Else If (noncolin) Then
            !
                nspin = 4
                !
            Else
            !
                nspin = 1
                !
            End If
            !
        End If
        !
        Call mp_bcast (lsda, ionode_id, intra_image_comm)
        Call mp_bcast (nspin, ionode_id, intra_image_comm)
        Call mp_bcast (noncolin, ionode_id, intra_image_comm)
        Call mp_bcast (npol, ionode_id, intra_image_comm)
        Call mp_bcast (lspinorb, ionode_id, intra_image_comm)
        Call mp_bcast (domag, ionode_id, intra_image_comm)
        !
        lspin_read = .True.
        !
        Return
        !
    End Subroutine read_spin
    !
    !--------------------------------------------------------------------------
    Subroutine read_magnetization (ierr)
    !------------------------------------------------------------------------
    !
        Use klist, Only: two_fermi_energies, nelup, neldw
        Use ener, Only: ef_up, ef_dw
        Use lsda_mod, Only: starting_magnetization
        Use noncollin_module, Only: angle1, angle2, i_cons, mcons, bfield, lambda
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        Logical :: found
        Integer :: i, nsp
        !
        ierr = 0
        If (lstarting_mag_read) Return
        !
        !
        If (ionode) Then
        !
            Call qexml_read_magnetization (starting_magnetization=starting_magnetization, angle1=angle1, angle2=angle2, two_fermi_energies=two_fermi_energies, i_cons=i_cons, mcons=mcons, bfield=bfield, ef_up=ef_up, ef_dw=ef_dw, nelup=nelup, neldw=neldw, lambda=lambda, found=found, ierr=ierr)
            !
            angle1 (:) = angle1 (:) * PI / 180.d0
            angle2 (:) = angle2 (:) * PI / 180.d0
            !
            If (two_fermi_energies) Then
            !
                ef_up = ef_up * e2
                ef_dw = ef_dw * e2
                !
            End If
            !
        End If
        !
        Call mp_bcast (found, ionode_id, intra_image_comm)
        !
        If (found) Then
        !
            Call mp_bcast (starting_magnetization, ionode_id, intra_image_comm)
            Call mp_bcast (angle1, ionode_id, intra_image_comm)
            Call mp_bcast (angle2, ionode_id, intra_image_comm)
            Call mp_bcast (two_fermi_energies, ionode_id, intra_image_comm)
            Call mp_bcast (i_cons, ionode_id, intra_image_comm)
            Call mp_bcast (mcons, ionode_id, intra_image_comm)
            Call mp_bcast (bfield, ionode_id, intra_image_comm)
            Call mp_bcast (nelup, ionode_id, intra_image_comm)
            Call mp_bcast (neldw, ionode_id, intra_image_comm)
            Call mp_bcast (ef_up, ionode_id, intra_image_comm)
            Call mp_bcast (ef_dw, ionode_id, intra_image_comm)
            Call mp_bcast (lambda, ionode_id, intra_image_comm)
            !
        End If
        !
        lstarting_mag_read = .True.
        !
        Return
        !
    End Subroutine read_magnetization
    !
    !------------------------------------------------------------------------
    Subroutine read_xc (ierr)
    !------------------------------------------------------------------------
    !
        Use ions_base, Only: nsp
        Use funct, Only: enforce_input_dft
        Use ldaU, Only: lda_plus_u, lda_plus_u_kind, Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_J, Hubbard_alpha, Hubbard_J0, Hubbard_beta, U_projection
        Use kernel_table, Only: vdw_table_name
        Use acfdt_ener, Only: acfdt_in_pw
        Use control_flags, Only: llondon, lxdm, ts_vdw
        Use london_module, Only: scal6, lon_rcut
        Use tsvdw_module, Only: vdw_isolated
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        !
        Character (Len=20) :: dft_name
        Integer :: nsp_, inlc
        Logical :: nomsg = .True.
        !
        ierr = 0
        If (lxc_read) Return
        !
        If ( .Not. lions_read) Call errore ('read_xc', 'read ions first', 1)
        !
        If (ionode) Then
        !
            Call qexml_read_xc (dft_name, lda_plus_u, lda_plus_u_kind, U_projection, Hubbard_lmax, Hubbard_l, nsp_, Hubbard_U, Hubbard_J, Hubbard_J0, Hubbard_alpha, Hubbard_beta, inlc, vdw_table_name, acfdt_in_pw, llondon, scal6, lon_rcut, lxdm, ts_vdw, vdw_isolated, ierr)
            !
        End If
        !
        Call mp_bcast (dft_name, ionode_id, intra_image_comm)
        Call mp_bcast (lda_plus_u, ionode_id, intra_image_comm)
        Call mp_bcast (inlc, ionode_id, intra_image_comm)
        Call mp_bcast (llondon, ionode_id, intra_image_comm)
        Call mp_bcast (lxdm, ionode_id, intra_image_comm)
        Call mp_bcast (ts_vdw, ionode_id, intra_image_comm)
        !
        If (lda_plus_u) Then
        !
            Call mp_bcast (lda_plus_u_kind, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_lmax, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_l, ionode_id, intra_image_comm)
            Call mp_bcast (U_projection, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_U, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_J, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_J0, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_alpha, ionode_id, intra_image_comm)
            Call mp_bcast (Hubbard_beta, ionode_id, intra_image_comm)
            !
        End If
        !
        If (llondon) Then
            Call mp_bcast (scal6, ionode_id, intra_image_comm)
            Call mp_bcast (lon_rcut, ionode_id, intra_image_comm)
        End If
        !
        If (ts_vdw) Then
            Call mp_bcast (vdw_isolated, ionode_id, intra_image_comm)
        End If
        !
        ! SCF EXX/RPA
        !
        Call mp_bcast (acfdt_in_pw, ionode_id, intra_image_comm)
        !
        If (acfdt_in_pw) dft_name = 'NOX-NOC'
        !
        If (inlc > 0) Then
            Call mp_bcast (vdw_table_name, ionode_id, intra_image_comm)
        End If
        !
        If (llondon) Then
            Call mp_bcast (scal6, ionode_id, intra_image_comm)
            Call mp_bcast (lon_rcut, ionode_id, intra_image_comm)
        End If
        !
        If (ts_vdw) Then
            Call mp_bcast (vdw_isolated, ionode_id, intra_image_comm)
        End If
        !
        ! SCF EXX/RPA
        !
        Call mp_bcast (acfdt_in_pw, ionode_id, intra_image_comm)
        !
        If (acfdt_in_pw) dft_name = 'NOX-NOC'
        !
        ! discard any further attempt to set a different dft
        Call enforce_input_dft (dft_name, nomsg)
        !
        lxc_read = .True.
        !
        Return
        !
    End Subroutine read_xc
    !
    !
    !------------------------------------------------------------------------
    Subroutine read_brillouin_zone (ierr)
    !------------------------------------------------------------------------
    !
        Use lsda_mod, Only: lsda
        Use klist, Only: nkstot, xk, wk, qnorm
        Use start_k, Only: nks_start, xk_start, wk_start, nk1, nk2, nk3, k1, k2, k3
        Use symm_base, Only: nrot, s, sname
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        Character (iotk_attlenx) :: attr
        !
        Integer :: i, ik, num_k_points
        Logical :: found
        !
        ierr = 0
        If (lbz_read) Return
        !
        !
        If (ionode) Then
        !
        ! xk_start and wk_start are ALLOCATABLE inside the function
            Call qexml_read_bz (num_k_points=num_k_points, xk=xk, wk=wk, k1=k1, k2=k2, k3=k3, nk1=nk1, nk2=nk2, nk3=nk3, nks_start=nks_start, xk_start=xk_start, wk_start=wk_start, qnorm=qnorm, ierr=ierr)
            !
            nkstot = num_k_points
            !
            If (lsda) nkstot = num_k_points * 2
            !
            Do ik = 1, num_k_points
            !
                If (lsda) Then
                !
                    xk (:, ik+num_k_points) = xk (:, ik)
                    !
                    wk (ik+num_k_points) = wk (ik)
                    !
                End If
                !
            End Do
            !
        End If
        !
        Call mp_bcast (nkstot, ionode_id, intra_image_comm)
        Call mp_bcast (xk, ionode_id, intra_image_comm)
        Call mp_bcast (wk, ionode_id, intra_image_comm)
        Call mp_bcast (nk1, ionode_id, intra_image_comm)
        Call mp_bcast (nk2, ionode_id, intra_image_comm)
        Call mp_bcast (nk3, ionode_id, intra_image_comm)
        Call mp_bcast (k1, ionode_id, intra_image_comm)
        Call mp_bcast (k2, ionode_id, intra_image_comm)
        Call mp_bcast (k3, ionode_id, intra_image_comm)
        Call mp_bcast (qnorm, ionode_id, intra_image_comm)
        !
        Call mp_bcast (nks_start, ionode_id, intra_image_comm)
        If (nks_start > 0 .And. .Not. ionode) Then
            If ( .Not. ALLOCATED(xk_start)) ALLOCATE (xk_start(3, nks_start))
            If ( .Not. ALLOCATED(wk_start)) ALLOCATE (wk_start(nks_start))
        End If
        If (nks_start > 0) Then
            Call mp_bcast (xk_start, ionode_id, intra_image_comm)
            Call mp_bcast (wk_start, ionode_id, intra_image_comm)
        End If
        Call mp_bcast (nrot, ionode_id, intra_image_comm)
        Call mp_bcast (s, ionode_id, intra_image_comm)
        Call mp_bcast (sname, ionode_id, intra_image_comm)
        !
        lbz_read = .True.
        !
        Return
        !
    End Subroutine read_brillouin_zone
    !
    !------------------------------------------------------------------------
    Subroutine read_occupations (ierr)
    !------------------------------------------------------------------------
    !
        Use lsda_mod, Only: lsda, nspin
        Use fixed_occ, Only: tfixed_occ, f_inp
        Use ktetra, Only: ntetra, tetra
        Use klist, Only: ltetra, lgauss, ngauss, degauss, smearing
        Use electrons_base, Only: nupdwn
        Use wvfct, Only: nbnd
        !
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        Character (iotk_attlenx) :: attr
        !
        Integer :: i
        Logical :: found
        !
        ierr = 0
        If (locc_read) Return
        !
        If (ionode) Then
        !
        ! necessary to don't send nbnd and nspin as input in read_occ
            If ( .Not. ALLOCATED(f_inp)) Then
            !
                If (nspin == 4) Then
                    Allocate (f_inp(nbnd, 1))
                Else
                    Allocate (f_inp(nbnd, nspin))
                End If
                !
            End If
            !
            f_inp (:, :) = 0.0d0
            !
            Call qexml_read_occ (lgauss=lgauss, ngauss=ngauss, degauss=degauss, ltetra=ltetra, ntetra=ntetra, tetra=tetra, tfixed_occ=tfixed_occ, NSTATES_UP=nupdwn(1), NSTATES_DW=nupdwn(2), INPUT_OCC=f_inp, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            If (lgauss) Then
            !
                Select Case (ngauss)
                Case (0)
                    smearing = 'gaussian'
                Case (1)
                    smearing = 'Methfessel-Paxton'
                Case (-1)
                    smearing = 'Marzari-Vanderbilt'
                Case (-99)
                    smearing = 'Fermi-Dirac'
                Case Default
                    Call errore ('read_occupations', 'wrong smearing index', Abs(1000+ngauss))
                End Select
                !
                degauss = degauss * e2
                !
            Else
            !
                ngauss = 0
                degauss = 0.d0
                !
            End If
            !
            If ( .Not. ltetra) Then
            !
                ntetra = 0
                !
            End If
            !
            If ( .Not. tfixed_occ) Then
            !
                Deallocate (f_inp)
                !
            End If
            !
            !
        End If
        !
        Call mp_bcast (lgauss, ionode_id, intra_image_comm)
        !
        If (lgauss) Then
        !
            Call mp_bcast (ngauss, ionode_id, intra_image_comm)
            Call mp_bcast (degauss, ionode_id, intra_image_comm)
            Call mp_bcast (smearing, ionode_id, intra_image_comm)
            !
        End If
        !
        Call mp_bcast (ltetra, ionode_id, intra_image_comm)
        !
        If (ltetra) Then
        !
            Call mp_bcast (ntetra, ionode_id, intra_image_comm)
            Call mp_bcast (tetra, ionode_id, intra_image_comm)
            !
        End If
        !
        Call mp_bcast (tfixed_occ, ionode_id, intra_image_comm)
        !
        If (tfixed_occ) Then
        !
            Call mp_bcast (nupdwn, ionode_id, intra_image_comm)
            !
            If ( .Not. ALLOCATED(f_inp)) Then
            !
                If (nspin == 4) Then
                    Allocate (f_inp(nbnd, 1))
                Else
                    Allocate (f_inp(nbnd, nspin))
                End If
                !
            End If
            !
            Call mp_bcast (f_inp, ionode_id, intra_image_comm)
            !
        End If
        !
        locc_read = .True.
        !
        Return
        !
    End Subroutine read_occupations
    !
    !
    !
    !------------------------------------------------------------------------
    Subroutine read_band_structure (dirname, ierr)
    !------------------------------------------------------------------------
    !
        Use control_flags, Only: lkpoint_dir
        Use basis, Only: natomwfc
        Use lsda_mod, Only: lsda, isk
        Use klist, Only: nkstot, wk, nelec
        Use wvfct, Only: et, wg, nbnd
        Use ener, Only: ef, ef_up, ef_dw
        !
        Implicit None
        !
        Character (Len=*), Intent (In) :: dirname
        Integer, Intent (Out) :: ierr
        !
        Integer :: ik, ik_eff, num_k_points
        Logical :: found, two_fermi_energies_
        Character (Len=256) :: filename
        !
        ierr = 0
        If (lbs_read) Return
        !
        If ( .Not. lspin_read) Call errore ('read_band_structure', 'read spin first', 1)
        If ( .Not. lbz_read) Call errore ('read_band_structure', 'read band_structure first', 1)
        !
        !
        If (ionode) Then
        ! we don't need to read nspin, noncolin
            Call qexml_read_bands_info (nbnd=nbnd, num_k_points=num_k_points, natomwfc=natomwfc, nelec=nelec, ef=ef, two_fermi_energies=two_fermi_energies_, ef_up=ef_up, ef_dw=ef_dw, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
            If ( .Not. two_fermi_energies_) Then
                ef = ef * e2
            Else
                ef = 0.d0
                ef_up = ef_up * e2
                ef_dw = ef_dw * e2
            End If
            !
        End If
        !
        num_k_points = nkstot
        !
        If (lsda) num_k_points = nkstot / 2
        !
        If (ionode) Then
        !
            If ( .Not. lkpoint_dir) filename = TRIM (dirname) // '/' // TRIM (xmlpun) // '.eig'
            !
            Call qexml_read_bands_pw (num_k_points, nbnd, nkstot, lsda, lkpoint_dir, filename, isk=isk, et=et, wg=wg, ierr=ierr)
            !
            et (:, :) = et (:, :) * e2
            !
            Forall (ik=1:nkstot) wg (:, ik) = wg (:, ik) * wk (ik)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        Call mp_bcast (nelec, ionode_id, intra_image_comm)
        Call mp_bcast (natomwfc, ionode_id, intra_image_comm)
        Call mp_bcast (nbnd, ionode_id, intra_image_comm)
        Call mp_bcast (isk, ionode_id, intra_image_comm)
        Call mp_bcast (et, ionode_id, intra_image_comm)
        Call mp_bcast (wg, ionode_id, intra_image_comm)
        Call mp_bcast (ef, ionode_id, intra_image_comm)
        !
        lbs_read = .True.
        !
        Return
        !
    End Subroutine read_band_structure
    !
    !------------------------------------------------------------------------
    Subroutine read_wavefunctions (dirname, ierr)
    !------------------------------------------------------------------------
    !
    ! ... This routines reads wavefunctions from the new file format and
    ! ... writes them into the old format
    !
        Use control_flags, Only: twfcollect, lkpoint_dir
        Use cell_base, Only: tpiba2
        Use lsda_mod, Only: nspin, isk
        Use klist, Only: nkstot, wk, nks, xk, ngk
        Use wvfct, Only: npw, npwx, et, wg, nbnd
        Use gvecw, Only: ecutwfc
        Use wavefunctions_module, Only: evc
        Use io_files, Only: nwordwfc, iunwfc
        Use buffers, Only: save_buffer
        Use gvect, Only: ngm, ngm_g, g, ig_l2g
        Use noncollin_module, Only: noncolin, npol
        Use mp_images, Only: nproc_image, intra_image_comm
        Use mp_pools, Only: kunit, nproc_pool, me_pool, root_pool, intra_pool_comm, inter_pool_comm
        Use mp_bands, Only: me_bgrp, nbgrp, root_bgrp, intra_bgrp_comm
        !
#if defined __HDF5
        Use hdf5_qe, Only: evc_hdf5_write, read_attributes_hdf5, prepare_for_reading_final
        Use mp_pools, Only: inter_pool_comm
#endif
!
        Implicit None
        !
        Character (Len=*), Intent (In) :: dirname
        Integer, Intent (Out) :: ierr
        !
        Character (Len=256) :: filename
        Integer :: ik, ipol, ik_eff, num_k_points
        Integer, Allocatable :: kisort (:)
        Integer :: nkl, nkr, npwx_g
        Integer :: nupdwn (2), ike, iks, npw_g, ispin
        Integer, External :: global_kpoint_index
        Integer, Allocatable :: ngk_g (:)
        Integer, Allocatable :: igk_l2g (:, :), igk_l2g_kdip (:, :)
        Logical :: opnd
        Real (DP), Allocatable :: gk (:)
        Real (DP) :: scalef
        !
        !
        ! The ierr output var is actually not given any value
        ! except this initialization
        !
        ierr = 0
        !
        iks = global_kpoint_index (nkstot, 1)
        ike = iks + nks - 1
        !
        ! ... find out the global number of G vectors: ngm_g
        !
        ngm_g = ngm
        !
        Call mp_sum (ngm_g, intra_bgrp_comm)
        !
        ! ... build the igk_l2g array, yielding the correspondence between
        ! ... the local k+G index and the global G index - see also ig_l2g
        !
        Allocate (igk_l2g(npwx, nks))
        igk_l2g = 0
        !
        Allocate (kisort(npwx), gk(npwx))
        !
        Do ik = 1, nks
        !
            kisort = 0
            npw = npwx
            !
            Call gk_sort (xk(1, ik+iks-1), ngm, g, ecutwfc/tpiba2, npw, kisort(1), gk)
            !
            Call gk_l2gmap (ngm, ig_l2g(1), npw, kisort(1), igk_l2g(1, ik))
            !
            ngk (ik) = npw
            !
        End Do
        !
        Deallocate (gk, kisort)
        !
        ! ... compute the global number of G+k vectors for each k point
        !
        Allocate (ngk_g(nkstot))
        !
        ngk_g = 0
        ngk_g (iks:ike) = ngk (1:nks)
        !
        Call mp_sum (ngk_g, inter_pool_comm)
        Call mp_sum (ngk_g, intra_pool_comm)
        ngk_g = ngk_g / nbgrp
        !
        ! ... compute the Maximum G vector index among all G+k an processors
        !
        npw_g = MAXVAL (igk_l2g(:, :))
        !
        Call mp_max (npw_g, inter_pool_comm)
        Call mp_max (npw_g, intra_pool_comm)
        !
        !
        ! ... compute the Maximum number of G vector among all k points
        !
        npwx_g = MAXVAL (ngk_g(1:nkstot))
        !
        !
        ! ... define a further l2g map to read gkvectors and wfc coherently
        !
        Allocate (igk_l2g_kdip(npwx_g, nks))
        igk_l2g_kdip = 0
        !
        Do ik = iks, ike
        !
            Call gk_l2gmap_kdip (npw_g, ngk_g(ik), ngk(ik-iks+1), igk_l2g(1, ik-iks+1), igk_l2g_kdip(1, ik-iks+1))
        End Do
        !
        !
        If (ionode) Then
        !
            Call iotk_scan_begin (iunpun, "EIGENVECTORS")
            !
        End If
        !
        num_k_points = nkstot
        !
        If (nspin == 2) num_k_points = nkstot / 2
        !
        k_points_loop: Do ik = 1, num_k_points
        !
            If (ionode) Then
            !
                Call iotk_scan_begin (iunpun, "K-POINT"// TRIM(iotk_index(ik)))
                !
                If (nspin == 2 .Or. noncolin) Then
                !
                    Call iotk_scan_begin (iunpun, "WFC.1", found=twfcollect)
                    If (twfcollect) Call iotk_scan_end (iunpun, "WFC.1")
                    !
                Else
                !
                    Call iotk_scan_begin (iunpun, "WFC", found=twfcollect)
                    If (twfcollect) Call iotk_scan_end (iunpun, "WFC")
                    !
                End If
                !
            End If
            !
            Call mp_bcast (twfcollect, ionode_id, intra_image_comm)
            !
            If ( .Not. twfcollect) Then
            !
                If (ionode) Then
                !
                    Call iotk_scan_end (iunpun, "K-POINT"// TRIM(iotk_index(ik)))
                    !
                End If
                !
                Exit k_points_loop
                !
            End If
            !
            If (nspin == 2) Then
            !
                ispin = 1
                evc = (0.0_DP, 0.0_DP)
                !
                ! ... no need to read isk here: they are read from band structure
                ! ... and correctly distributed across pools in read_file
                !!! isk(ik) = 1
                !
                If (ionode) Then
                !
                    filename = TRIM (qexml_wfc_filename(dirname, 'evc', ik, ispin, DIR=lkpoint_dir))
                    !
                End If
                !
                Call read_wfc (iunout, ik, nkstot, kunit, ispin, nspin, evc, npw_g, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, scalef, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                !
                If ((ik >= iks) .And. (ik <= ike)) Then
                !
                    Call save_buffer (evc, nwordwfc, iunwfc, (ik-iks+1))
                    !
                End If
                !
                ispin = 2
                ik_eff = ik + num_k_points
                evc = (0.0_DP, 0.0_DP)
                !
                ! ... no need to read isk here (see above why)
                !isk(ik_eff) = 2
                !
                If (ionode) Then
                !
                    filename = TRIM (qexml_wfc_filename(dirname, 'evc', ik, ispin, DIR=lkpoint_dir))
                    !
                End If
                !
                Call read_wfc (iunout, ik_eff, nkstot, kunit, ispin, nspin, evc, npw_g, nbnd, igk_l2g_kdip(:, ik_eff-iks+1), ngk(ik_eff-iks+1), filename, scalef, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                !
                If ((ik_eff >= iks) .And. (ik_eff <= ike)) Then
                !
                    Call save_buffer (evc, nwordwfc, iunwfc, (ik_eff-iks+1))
                    !
                End If
                !
            Else
            !
            ! ... no need to read isk here (see above why)
            !isk(ik) = 1
            !
                evc = (0.0_DP, 0.0_DP)
                If (noncolin) Then
                !
                    Do ipol = 1, npol
                    !
                        If (ionode) Then
                        !
                            filename = TRIM (qexml_wfc_filename(dirname, 'evc', ik, ipol, DIR=lkpoint_dir))
                            !
                        End If
                        !
                        !!! TEMP
                        nkl = (ipol-1) * npwx + 1
                        nkr = ipol * npwx
                        Call read_wfc (iunout, ik, nkstot, kunit, ispin, npol, evc(nkl:nkr, :), npw_g, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, scalef, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                        !
                    End Do
                    !
                Else
                !
                    If (ionode) Then
                    !
                        filename = TRIM (qexml_wfc_filename(dirname, 'evc',  ik, DIR=lkpoint_dir))
                        !
                    End If
                    !
                    ! workaround for pot parallelization ( Viet Nguyen / SdG )
                    ! -pot parallelization uses mp_image communicators
                    ! note that ionode must be also reset in the similar way
                    ! to image parallelization
                    Call read_wfc (iunout, ik, nkstot, kunit, ispin, nspin, evc, npw_g, nbnd, igk_l2g_kdip(:, ik-iks+1), ngk(ik-iks+1), filename, scalef, ionode, root_pool, intra_pool_comm, inter_pool_comm, intra_image_comm)
                    !
                End If
                !
                If ((ik >= iks) .And. (ik <= ike)) Then
                !
                    Call save_buffer (evc, nwordwfc, iunwfc, (ik-iks+1))
                    !
                    ! the following two line can be used to debug read_wfc
                    ! WRITE(200+10*ik+me_pool,fmt="(2D18.10)") evc
                    ! CLOSE(200+10*ik+me_pool )
                    !
                End If
                !
            End If
            !
            If (ionode) Then
            !
                Call iotk_scan_end (iunpun, "K-POINT"// TRIM(iotk_index(ik)))
                !
            End If
            !
        End Do k_points_loop
        !
        Deallocate (igk_l2g)
        Deallocate (igk_l2g_kdip)
        !
        If (ionode) Then
        !
            Call iotk_scan_end (iunpun, "EIGENVECTORS")
            !
            !CALL iotk_close_read( iunpun )
            !
        End If
        !
        Return
        !
    End Subroutine read_wavefunctions
    !
    !------------------------------------------------------------------------
    Subroutine read_ef (ierr)
    !------------------------------------------------------------------------
    !
    ! ... this routine reads the Fermi energy and the number of electrons
    !
        Use ener, Only: ef, ef_up, ef_dw
        Use klist, Only: two_fermi_energies, nelec
        !
        Implicit None
        Integer, Intent (Out) :: ierr
        !
        ! ... then selected tags are read from the other sections
        !
        If (ionode) Then
        !
            Call qexml_read_bands_info (ef=ef, ef_up=ef_up, ef_dw=ef_dw, two_fermi_energies=two_fermi_energies, nelec=nelec, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            If ( .Not. two_fermi_energies) Then
                ef = ef * e2
                ef_up = 0.d0
                ef_dw = 0.d0
            Else
                ef = 0.d0
                ef_up = ef_up * e2
                ef_dw = ef_dw * e2
            End If
            !
        End If
        !
        Call mp_bcast (two_fermi_energies, ionode_id, intra_image_comm)
        Call mp_bcast (ef, ionode_id, intra_image_comm)
        Call mp_bcast (ef_up, ionode_id, intra_image_comm)
        Call mp_bcast (ef_dw, ionode_id, intra_image_comm)
        Call mp_bcast (nelec, ionode_id, intra_image_comm)
        !
        Return
        !
    End Subroutine read_ef
    !
    !------------------------------------------------------------------------
    Subroutine read_exx (ierr)
    !------------------------------------------------------------------------
    !
    ! ... read EXX variables
    !
        Use funct, Only: set_exx_fraction, set_screening_parameter, set_gau_parameter, enforce_input_dft, start_exx
        Use exx, Only: x_gamma_extrapolation, nq1, nq2, nq3, exxdiv_treatment, yukawa, ecutvcut, ecutfock
        Implicit None
        !
        Integer, Intent (Out) :: ierr
        Real (DP) :: exx_fraction, screening_parameter, gau_parameter
        Logical :: exx_is_active, found
        !
        If (ionode) Then
            Call qexml_read_exx (x_gamma_extrapolation=x_gamma_extrapolation, NQX1=nq1, NQX2=nq2, NQX3=nq3, exxdiv_treatment=exxdiv_treatment, yukawa=yukawa, ecutvcut=ecutvcut, exx_fraction=exx_fraction, screening_parameter=screening_parameter, gau_parameter=gau_parameter, exx_is_active=exx_is_active, ecutfock=ecutfock, found=found, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Return
        !
        Call mp_bcast (found, ionode_id, intra_image_comm)
        !
        If ( .Not. found) Return
        !
        Call mp_bcast (x_gamma_extrapolation, ionode_id, intra_image_comm)
        Call mp_bcast (nq1, ionode_id, intra_image_comm)
        Call mp_bcast (nq2, ionode_id, intra_image_comm)
        Call mp_bcast (nq3, ionode_id, intra_image_comm)
        Call mp_bcast (exxdiv_treatment, ionode_id, intra_image_comm)
        Call mp_bcast (yukawa, ionode_id, intra_image_comm)
        Call mp_bcast (ecutvcut, ionode_id, intra_image_comm)
        Call mp_bcast (exx_fraction, ionode_id, intra_image_comm)
        Call mp_bcast (screening_parameter, ionode_id, intra_image_comm)
        Call mp_bcast (gau_parameter, ionode_id, intra_image_comm)
        Call mp_bcast (exx_is_active, ionode_id, intra_image_comm)
        Call mp_bcast (ecutfock, ionode_id, intra_image_comm)
        !
        Call set_exx_fraction (exx_fraction)
        Call set_screening_parameter (screening_parameter)
        Call set_gau_parameter (gau_parameter)
        If (exx_is_active) Call start_exx ()
        !
        Return
        !
    End Subroutine read_exx
    !
    !------------------------------------------------------------------------
    Subroutine read_esm (ierr)
    !------------------------------------------------------------------------
    !
    ! ... this routine reads only nelec and ef
    !
        Use esm, Only: esm_nfit, esm_efield, esm_w, esm_a, esm_bc
        !
        Implicit None
        Integer, Intent (Out) :: ierr
        !
        ! ... then selected tags are read from the other sections
        !
        If (ionode) Then
        !
            Call qexml_read_esm (esm_nfit=esm_nfit, esm_efield=esm_efield, esm_w=esm_w, esm_a=esm_a, esm_bc=esm_bc, ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        If (ierr > 0) Return
        !
        Call mp_bcast (esm_nfit, ionode_id, intra_image_comm)
        Call mp_bcast (esm_efield, ionode_id, intra_image_comm)
        Call mp_bcast (esm_w, ionode_id, intra_image_comm)
        Call mp_bcast (esm_a, ionode_id, intra_image_comm)
        Call mp_bcast (esm_bc, ionode_id, intra_image_comm)
        !
        Return
        !
    End Subroutine read_esm
    !
    !------------------------------------------------------------------------
    Subroutine read_ (dirname, ierr)
    !------------------------------------------------------------------------
    !
    ! ... this is a template for a "read section" subroutine
    !
        Implicit None
        !
        Character (Len=*), Intent (In) :: dirname
        Integer, Intent (Out) :: ierr
        !
        Integer :: idum
        !
        !
        If (ionode) Then
        !
            Call iotk_open_read (iunpun, FILE=TRIM(dirname)//'/'// TRIM(xmlpun), ierr=ierr)
            !
        End If
        !
        Call mp_bcast (ierr, ionode_id, intra_image_comm)
        !
        If (ierr > 0) Return
        !
        If (ionode) Then
        !
            Call iotk_scan_begin (iunpun, "")
            !
            Call iotk_scan_end (iunpun, "")
            !
            Call iotk_close_read (iunpun)
            !
        End If
        !
        Call mp_bcast (idum, ionode_id, intra_image_comm)
        !
        Return
        !
    End Subroutine read_
    !
    !----------------------------------------------------------------------------
    Subroutine gk_l2gmap (ngm, ig_l2g, ngk, igk, igk_l2g)
    !----------------------------------------------------------------------------
    !
    ! ... This subroutine maps local G+k index to the global G vector index
    ! ... the mapping is used to collect wavefunctions subsets distributed
    ! ... across processors.
    ! ... Written by Carlo Cavazzoni
    !
        Implicit None
        !
        ! ... Here the dummy variables
        !
        Integer, Intent (In) :: ngm, ngk, igk (ngk), ig_l2g (ngm)
        Integer, Intent (Out) :: igk_l2g (ngk)
        Integer :: ig
        !
        ! ... input: mapping between local and global G vector index
        !
        Do ig = 1, ngk
        !
            igk_l2g (ig) = ig_l2g (igk(ig))
            !
        End Do
        !
        Return
        !
    End Subroutine gk_l2gmap
    !
    !-----------------------------------------------------------------------
    Subroutine gk_l2gmap_kdip (npw_g, ngk_g, ngk, igk_l2g, igk_l2g_kdip, igwk)
    !-----------------------------------------------------------------------
    !
    ! ... This subroutine maps local G+k index to the global G vector index
    ! ... the mapping is used to collect wavefunctions subsets distributed
    ! ... across processors.
    ! ... This map is used to obtained the G+k grids related to each kpt
    !
        Implicit None
        !
        ! ... Here the dummy variables
        !
        Integer, Intent (In) :: npw_g, ngk_g, ngk
        Integer, Intent (In) :: igk_l2g (ngk)
        Integer, Optional, Intent (Out) :: igwk (ngk_g), igk_l2g_kdip (ngk)
        !
        Integer, Allocatable :: igwk_ (:), itmp (:), igwk_lup (:)
        Integer :: ig, ig_, ngg
        !
        !
        Allocate (itmp(npw_g))
        Allocate (igwk_(ngk_g))
        !
        itmp (:) = 0
        igwk_ (:) = 0
        !
        !
        Do ig = 1, ngk
        !
            itmp (igk_l2g(ig)) = igk_l2g (ig)
            !
        End Do
        !
        !
        Call mp_sum (itmp, intra_bgrp_comm)
        !
        ngg = 0
        Do ig = 1, npw_g
        !
            If (itmp(ig) == ig) Then
            !
                ngg = ngg + 1
                !
                igwk_ (ngg) = ig
                !
            End If
            !
        End Do
        !
        If (ngg /= ngk_g) Call errore ('gk_l2gmap_kdip', 'unexpected dimension in ngg', 1)
        !
        If (PRESENT(igwk)) Then
        !
            igwk (1:ngk_g) = igwk_ (1:ngk_g)
            !
        End If
        !
        If (PRESENT(igk_l2g_kdip)) Then
        !
            Allocate (igwk_lup(npw_g))
            !
            !$omp parallel private(ig_, ig)
            !$omp workshare
            igwk_lup = 0
            !$omp end workshare
            !$omp do
            Do ig_ = 1, ngk_g
                igwk_lup (igwk_(ig_)) = ig_
            End Do
            !$omp end do
            !$omp do
            Do ig = 1, ngk
                igk_l2g_kdip (ig) = igwk_lup (igk_l2g(ig))
            End Do
            !$omp end do
            !$omp end parallel
            !
            Deallocate (igwk_lup)
            !
        End If
        !
        Deallocate (itmp, igwk_)
        !
        Return
        !
    End Subroutine gk_l2gmap_kdip
#endif
!
End Module pw_restart
!