wannier90_interface.F90

Go to the documentation of this file.

!! Copyright (C) 2017-2019 H. Huebener, N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
program wannier90_interface
  use batch_oct_m
  use calc_mode_par_oct_m
  use comm_oct_m
  use command_line_oct_m
  use cube_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use fft_oct_m
  use global_oct_m
  use grid_oct_m
  use io_oct_m
  use io_binary_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use loct_oct_m
  use loct_math_oct_m
  use mesh_batch_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use orbitalset_oct_m
  use parser_oct_m
  use profiling_oct_m
  use restart_oct_m
  use states_elec_calc_oct_m
  use electrons_oct_m
  use space_oct_m
  use string_oct_m
  use states_elec_oct_m
  use states_elec_io_oct_m
  use states_elec_restart_oct_m
  use submesh_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use wfs_elec_oct_m
  use ylm_wannier_oct_m
 
  implicit none
 
  integer              :: w90_what, w90_mode
  integer(int64)       :: w90_what_default
 
  integer              :: ierr
  integer              :: dim, idim
  integer              :: ii, nik, iter, nst
 
  type(restart_t)      :: restart
  type(electrons_t), pointer :: sys
  logical              :: w90_spinors, scdm_proj, w90_scdm
  integer              :: w90_nntot, w90_num_bands, w90_num_kpts   ! w90 input parameters
  integer, allocatable :: w90_nnk_list(:,:)                        !
  character(len=80)    :: w90_prefix                               ! w90 input file prefix
  integer              :: w90_num_wann                             ! input paramter
  real(real64), allocatable   :: w90_proj_centers(:,:)                    ! projections centers
  integer, allocatable :: w90_proj_lmr(:,:)                        ! definitions for real valued Y_lm*R_r
  integer :: w90_nproj                                             ! number of such projections
  integer, allocatable :: w90_spin_proj_component(:)               ! up/down flag
  real(real64), allocatable   :: w90_spin_proj_axis(:,:)                  ! spin axis (not implemented)
  integer              :: w90_num_exclude
  logical, allocatable :: exclude_list(:)                          ! list of excluded bands
  integer, allocatable :: band_index(:)                            ! band index after exclusion
  logical              :: read_td_states
  integer              :: w90_spin_channel
  logical              :: w90_bloch_sums = .false.                 
 
  ! scdm variables
  integer, allocatable :: jpvt(:)
  complex(real64), allocatable   :: uk(:,:,:)                                ! SCDM-Wannier gauge matrices U(k)
  complex(real64), allocatable   :: psi(:,:)
  complex(real64), allocatable   :: chi(:,:), chi_diag(:,:),chi2(:,:)
  real(real64), allocatable   :: chi_eigenval(:), occ_temp(:)
  real(real64)         :: scdm_mu, scdm_sigma, smear,  kvec(3), factor(3)
  integer :: ist, jst, ik, idir, sender
 
  integer(int64) :: how
 
  call global_init()
  call parser_init()
 
  call messages_init()
  call io_init()
 
  call profiling_init(global_namespace)
 
  call fft_all_init(global_namespace)
  call unit_system_init(global_namespace)
 
  call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
  sys => electrons_t(global_namespace, int(option__calculationmode__dummy, int32))
  call sys%init_parallelization(mpi_world)
 
  !%Variable Wannier90Prefix
  !%Type string
  !%Default w90
  !%Section Utilities::oct-wannier90
  !%Description
  !% Prefix for wannier90 files
  !%End
  call parse_variable(global_namespace, 'Wannier90Prefix', 'w90', w90_prefix)
  if (w90_prefix == 'w90') then
    message(1) = "oct-wannier90: the prefix is set by default to w90"
    call messages_info(1)
  end if
 
  !%Variable Wannier90Mode
  !%Type integer
  !%Default 0
  !%Section Utilities::oct-wannier90
  !%Description
  !% Specifies which stage of the Wannier90 interface to use
  !%Option none 0
  !% Nothing is done.
  !%Option w90_setup 1
  !% Writes parts of the wannier90 input file <tt>w90_prefix.win</tt> corresponding to
  !% the octopus inp file. Importantly it generates the correct form of Monkhorst-Pack mesh
  !% written to the file w90_kpoints that has to be used in a gs calculation of Octopus by
  !% as <tt> include w90_kpoints </tt> instead of the <tt>%KpointsGrid</tt> block.
  !%Option w90_output 2
  !% Generates the relevant files for a wannier90 run, specified by the variable <tt>W90_interface_files</tt>.
  !% This needs files previously generated
  !% by <tt>wannier90.x -pp w90 </tt>
  !%Option w90_wannier 3
  !% Parse the output of wannier90 to generate the Wannier states on the real-space grid.
  !% The states will be written in the folder wannier. By default, the states are written as
  !% binary files, similar to the Kohn-Sham states.
  !%
  !% Not implemented for spinor states.
  !%End
  call parse_variable(global_namespace, 'Wannier90Mode', 0, w90_mode)
 
  if (w90_mode == 0) then
    message(1) = "Wannier90Mode must be set to a value different from 0."
    call messages_fatal(1)
  end if
 
  !%Variable Wannier90Files
  !%Type flag
  !%Default w90_mmn + w90_amn + w90_eig
  !%Section Utilities::oct-wannier90
  !%Description
  !% Specifies which files to generate.
  !% Example: <tt>w90_mmn + w90_unk</tt>
  !%Option w90_mmn bit(1)
  !% (see Wannier90 documentation)
  !%Option w90_unk bit(2)
  !% (see Wannier90 documentation)
  !%Option w90_amn bit(3)
  !% (see Wannier90 documentation)
  !%Option w90_eig bit(4)
  !% Eigenvalues. See Wannier90 documentation for more details.
  !%Option w90_spn bit(5)
  !% Spin. See Wannier90 documentation for more details.
  !%End
  w90_what_default = option__wannier90files__w90_mmn + option__wannier90files__w90_amn + option__wannier90files__w90_eig
  if (sys%st%d%ispin == spinors) w90_what_default = w90_what_default + option__wannier90files__w90_spn
  call parse_variable(global_namespace, 'Wannier90Files', w90_what_default, w90_what)
 
  !%Variable Wannier90UseTD
  !%Type logical
  !%Default no
  !%Section Utilities::oct-wannier90
  !%Description
  !% By default oct-wannier90 uses the ground-state states to compute the necessary information.
  !% By setting this variable to yes, oct-wannier90 will use the TD states instead.
  !%End
  call parse_variable(global_namespace, 'Wannier90UseTD', .false., read_td_states)
 
  !%Variable Wannier90UseSCDM
  !%Type logical
  !%Default no
  !%Section Utilities::oct-wannier90
  !%Description
  !% By default oct-wannier90 uses the projection method to generate the .amn file.
  !% By setting this variable to yes, oct-wannier90 will use SCDM method instead.
  !%End
  call parse_variable(global_namespace, 'Wannier90UseSCDM', .false., w90_scdm)
  if (w90_scdm) then
    !%Variable SCDMsigma
    !%Type float
    !%Default 0.2
    !%Section Utilities::oct-wannier90
    !%Description
    !% Broadening of SCDM smearing function.
    !%End
    call parse_variable(global_namespace, 'SCDMsigma', 0.2_real64, scdm_sigma)
 
    !%Variable SCDMmu
    !%Type float
    !%Section Utilities::oct-wannier90
    !%Description
    !% Energy range up to which states are considered for SCDM.
    !%End
    call parse_variable(global_namespace, 'SCDMmu', m_huge, scdm_mu)
  end if
 
  if (sys%kpoints%use_symmetries) then
    message(1) = 'oct-wannier90: k-points symmetries are not allowed'
    call messages_fatal(1)
  end if
  if (sys%kpoints%use_time_reversal) then
    message(1) = 'oct-wannier90: time-reversal symmetry is not allowed'
    call messages_fatal(1)
  end if
  if (sys%kpoints%reduced%nshifts > 1) then
    message(1) = 'oct-wannier90: Wannier90 does not allow for multiple shifts of the k-point grid'
    call messages_fatal(1)
  end if
 
  if (sys%st%d%ispin /= unpolarized) then
    call messages_experimental("oct-wannier90 with SpinComponnents /= unpolarized")
  end if
 
  w90_spinors = .false.
 
  !%Variable Wannier90SpinChannel
  !%Type integer
  !%Section Utilities::oct-wannier90
  !%Description
  !% Spin channel used for the Wannierization
  !%End
  call parse_variable(global_namespace, 'Wannier90SpinChannel', 1, w90_spin_channel)
 
  ! "Correct" rlattice/klattice for Wannier90 (3D periodicity assumed here).
  factor = m_one
  do idir = sys%space%periodic_dim+1, sys%space%dim
    factor(idir) = m_two * sys%gr%box%bounding_box_l(idir)
  end do
  call sys%ions%latt%scale(factor)
 
  ! create setup files
  select case (w90_mode)
  case (option__wannier90mode__w90_setup)
    call wannier90_setup(sys%ions, sys%kpoints, sys%space)
 
    ! load states and calculate interface files
  case (option__wannier90mode__w90_output)
    call wannier90_output()
 
  case (option__wannier90mode__w90_wannier)
    !%Variable Wannier90ComputeBlochSums
    !%Type logical
    !%Default no
    !%Section Utilities::oct-wannier90
    !%Description
    !% Only for Wannier90Mode = w90_wannier
    !% By setting this variable to yes, oct-wannier90 will also output the Bloch sums
    !% of the Wannier orbitals. This creates one file per k-point.
    !%End
    call parse_variable(global_namespace, 'Wannier90ComputeBlochSums', .false., w90_bloch_sums)
 
    call read_wannier90_files()
 
    ! normal interface run
    call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type = type_cmplx, skip=exclude_list)
    if (read_td_states) then
      call restart%init(global_namespace, restart_td, restart_type_load, &
        sys%mc, ierr, sys%gr)
    else
      call restart%init(global_namespace, restart_gs, restart_type_load, &
        sys%mc, ierr, sys%gr)
    end if
 
    if (ierr == 0) then
      call states_elec_look(restart, nik, dim, nst, ierr)
      if (sys%st%d%ispin == spin_polarized) then
        nik = nik / 2
      end if
      if (dim == sys%st%d%dim .and. nik == sys%kpoints%reduced%npoints .and. nst >= sys%st%nst) then
        call states_elec_load(restart, global_namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
          ierr, iter, label = ": wannier90", skip=exclude_list)
      else
        write(message(1),'(a)') 'Restart structure not commensurate.'
        call messages_fatal(1)
      end if
    end if
    call restart%end()
 
    call generate_wannier_states(sys%space, sys%gr, sys%ions, sys%st, sys%kpoints)
  case default
    message(1) = "Wannier90Mode is set to an unsupported value."
    call messages_fatal(1)
  end select
 
  safe_deallocate_a(exclude_list)
  safe_deallocate_a(band_index)
  safe_deallocate_a(w90_nnk_list)
  safe_deallocate_a(w90_proj_centers)
  safe_deallocate_a(w90_proj_lmr)
 
  safe_deallocate_p(sys)
  call fft_all_end()
  call io_end()
  call profiling_end(global_namespace)
  call messages_end()
  call parser_end()
  call global_end()
 
contains
 
  ! --------------------------------------------------------------------------
  subroutine wannier90_setup(ions, kpoints, space)
    type(ions_t),      intent(in) :: ions
    type(kpoints_t),   intent(in) :: kpoints
    class(space_t),    intent(in) :: space
 
    character(len=80) :: filename
    integer :: w90_win, ia, axis(3), npath
 
    push_sub(wannier90_setup)
 
    assert(space%dim == 3)
 
    ! open win file
    filename = trim(adjustl(w90_prefix)) //'.win'
    w90_win = io_open(trim(filename), global_namespace, action='write')
 
    write(w90_win,'(a)') '# this file has been created by the Octopus wannier90 utility'
    write(w90_win,'(a)') ' '
 
    ! write direct lattice vectors (in angstrom)
    write(w90_win,'(a)') 'begin unit_cell_cart'
    write(w90_win,'(a)') 'Ang'
    do idim = 1,3
      write(w90_win,'(f13.8,f13.8,f13.8)') units_from_atomic(unit_angstrom, ions%latt%rlattice(1:3,idim))
    end do
    write(w90_win,'(a)') 'end unit_cell_cart'
    write(w90_win,'(a)') ' '
 
    write(w90_win,'(a)') 'begin atoms_frac'
    do ia = 1, ions%natoms
      write(w90_win,'(a,2x,f13.8,f13.8,f13.8)') trim(ions%atom(ia)%label), ions%latt%cart_to_red(ions%pos(:, ia))
    end do
    write(w90_win,'(a)') 'end atoms_frac'
    write(w90_win,'(a)') ' '
 
    ! This is a default value. In practice, one should use projections
    write(w90_win,'(a)') 'use_bloch_phases = .true.'
    write(w90_win,'(a)') ' '
 
    write(w90_win,'(a10,i4)') 'num_bands ', sys%st%nst
    write(w90_win,'(a9,i4)') 'num_wann ', sys%st%nst
    write(w90_win,'(a)') ' '
 
    if (sys%st%d%ispin == spinors) then
      write(w90_win,'(a)') 'spinors = .true.'
    end if
    if (sys%st%d%ispin == spin_polarized) then
      if (w90_spin_channel == 1) then
        write(w90_win, '(a)') 'spin = up'
      else if (w90_spin_channel == 2) then
        write(w90_win, '(a)') 'spin = down'
      else
        message(1) = 'Wannier90SpinChannel value is invalid.'
        call messages_fatal(1)
      end if
    end if
 
    ! This is for convenience. This is needed for plotting the Wannier states, if requested.
    write(w90_win,'(a)')  'write_u_matrices = .true.'
    write(w90_win,'(a)')  'write_xyz = .true.'
    write(w90_win,'(a)') ' '
 
    if (kpoints%reduced%npoints == 1) then
      write(w90_win,'(a)')  'gamma_only = .true.'
      write(w90_win,'(a)') ' '
    else
      if (.not. parse_is_defined(global_namespace, 'KPointsGrid')) then
        message(1) = 'oct-wannier90: need Monkhorst-Pack grid. Please specify %KPointsGrid'
        call messages_fatal(1)
      end if
 
      ! In case the user used also a k-point path, we ignore it
      npath = kpoints%nkpt_in_path()
 
      axis(1:3) = kpoints%nik_axis(1:3)
      assert(product(kpoints%nik_axis(1:3)) == kpoints%reduced%npoints - npath)
 
      write(w90_win,'(a8,i4,i4,i4)')  'mp_grid =', axis(1:3)
      write(w90_win,'(a)') ' '
      write(w90_win,'(a)')  'begin kpoints '
      ! Put a minus sign here for the wrong convention in Octopus
 
      do ii = 1, kpoints%reduced%npoints-npath
        write(w90_win,'(f13.8,f13.8,f13.8)') - kpoints%reduced%red_point(1:3,ii)
      end do
      write(w90_win,'(a)')  'end kpoints '
    end if
 
    call io_close(w90_win)
 
    pop_sub(wannier90_setup)
 
  end subroutine wannier90_setup
 
  ! --------------------------------------------------------------------------
  subroutine wannier90_output()
    integer :: ik_real
 
    push_sub(wannier90_output)
 
    call read_wannier90_files()
 
    ! normal interface run
    call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type = type_cmplx, skip=exclude_list)
    if (read_td_states) then
      call restart%init(global_namespace, restart_td, restart_type_load, &
        sys%mc, ierr, sys%gr)
    else
      call restart%init(global_namespace, restart_gs, restart_type_load, &
        sys%mc, ierr, sys%gr)
    end if
 
    if (ierr == 0) then
      call states_elec_look(restart, nik, dim, nst, ierr)
      if (sys%st%d%ispin == spin_polarized) then
        nik = nik / 2
      end if
      if (dim == sys%st%d%dim .and. nik == sys%kpoints%reduced%npoints .and. nst >= sys%st%nst) then
        call states_elec_load(restart, global_namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
          ierr, iter, label = ": wannier90", skip=exclude_list)
      else
        write(message(1),'(a)') 'Restart structure not commensurate.'
        call messages_fatal(1)
      end if
    end if
    call restart%end()
 
    if (w90_scdm) then
      nik = w90_num_kpts
      safe_allocate(jpvt(1:sys%gr%np_global*sys%st%d%dim))
      safe_allocate(psi(1:sys%gr%np, 1:sys%st%d%dim))
      safe_allocate(occ_temp(1:w90_num_bands))
 
      ! smear the states at gamma
      do ist = 1, w90_num_bands
        occ_temp(ist)= sys%st%occ(ist, 1)
        sys%st%occ(ist, 1)=m_half*erfc((sys%st%eigenval(ist, 1)-scdm_mu) / scdm_sigma)
      end do
 
      call zstates_elec_rrqr_decomposition(sys%st, sys%namespace, sys%gr, w90_num_bands, .true., 1, jpvt)
 
      ! reset occupations at gamma
      do ist = 1, w90_num_bands
        sys%st%occ(ist, 1) = occ_temp(ist)
      end do
 
      safe_allocate(uk(1:w90_num_bands, 1:w90_num_bands, 1:nik))
 
      ! auxiliary arrays for scdm procedure
      safe_allocate(chi(1:w90_num_bands, 1:w90_num_bands))
      safe_allocate(chi_diag(1:w90_num_bands, 1:w90_num_bands))
      safe_allocate(chi2(1:w90_num_bands, 1:w90_num_bands))
      safe_allocate(chi_eigenval(1:w90_num_bands))
 
      chi(1:w90_num_bands, 1:w90_num_bands) = m_zero
 
      do ik = 1, nik
        kvec(:) = sys%kpoints%reduced%point(:, ik)
 
        if (sys%st%d%ispin == spin_polarized) then
          ik_real = (ik-1)*2 + w90_spin_channel
        else
          ik_real = ik
        end if
 
        do ist = 1, w90_num_bands
          sender = -1
          if (state_kpt_is_local(sys%st, ist, ik_real)) then
            call states_elec_get_state(sys%st, sys%gr, ist, ik_real, psi)
            sender = sys%st%d%kpt%mpi_grp%rank
          endif
          call sys%st%d%kpt%mpi_grp%allreduce_inplace(sender, 1, mpi_integer, mpi_max)
          assert(sender >= 0)
          call sys%st%d%kpt%mpi_grp%bcast(psi, sys%gr%np*sys%st%d%dim, mpi_double_complex, sender)
          smear=m_half * erfc((sys%st%eigenval(ist, ik_real) - scdm_mu) / scdm_sigma)
          ! NOTE: here check for domain parallelization
          do jst = 1, w90_num_bands
            chi(ist, jst) = smear * conjg(psi(jpvt(jst), 1)) &
              * exp(m_zi * dot_product(sys%gr%x(1:3, jpvt(jst)), kvec(1:3)))
          end do
        end do
 
        ! loewdin orhtogonalization of chi.chi
        ! this can also be done with SVD, which might be more stable!?
        chi_diag = matmul(conjg(transpose(chi)), chi)
        call lalg_eigensolve(w90_num_bands, chi_diag, chi_eigenval)
        chi2 = conjg(transpose(chi_diag))
 
        !we need the eigenvalues to be >0
        if (any(chi_eigenval(:) .lt. m_zero)) then
          message(1) = 'SCDM Wannierization failed because chi matrix is'
          message(2) = 'ill conditioned. Try increasingin scdm_sigma and/or'
          message(3) = 'change scdm_mu.'
          call messages_fatal(3)
        end if
 
        do ist = 1, w90_num_bands
          chi_eigenval(ist) = m_one / sqrt(chi_eigenval(ist))
          chi2(ist, 1:w90_num_bands) = chi_eigenval(ist) * chi2(ist, 1:w90_num_bands)
        end do
        ! the loewdin result would be: matmul(chi_diag,chi2)
        ! to get the wannier gauge U(k) we multiply this with the original chi
        uk(:,:,ik) = matmul(chi, matmul(chi_diag,chi2))
 
      end do
 
      safe_deallocate_a(chi)
      safe_deallocate_a(psi)
      safe_deallocate_a(chi_diag)
      safe_deallocate_a(chi2)
      safe_deallocate_a(chi_eigenval)
      safe_deallocate_a(jpvt)
      safe_deallocate_a(psi)
      safe_deallocate_a(occ_temp)
    end if
 
    ! ---- actual interface work ----------
    if (bitand(w90_what, option__wannier90files__w90_mmn) /= 0) then
      call create_wannier90_mmn(sys%gr, sys%st)
    end if
 
    if (bitand(w90_what, option__wannier90files__w90_unk) /= 0) then
      call write_unk(sys%space, sys%gr, sys%st, formatted=.false.)
    end if
 
    if (bitand(w90_what, option__wannier90files__w90_amn) /= 0) then
      call create_wannier90_amn(sys%space, sys%gr, sys%ions%latt, sys%st, sys%kpoints)
    end if
 
    if (bitand(w90_what, option__wannier90files__w90_eig) /= 0) then
      call create_wannier90_eig()
    end if
 
    if (bitand(w90_what, option__wannier90files__w90_spn) /= 0) then
      call create_wannier90_spn(sys%gr, sys%st)
    end if
 
    safe_deallocate_a(uk)
    safe_deallocate_a(w90_spin_proj_component)
    safe_deallocate_a(w90_spin_proj_axis)
 
    pop_sub(wannier90_output)
  end subroutine wannier90_output
 
  ! --------------------------------------------------------------------------
  subroutine read_wannier90_files()
    integer ::  w90_nnkp, itemp, dummyint, io, spin_channel_win
    character(len=80) :: filename, dummy, dummy1, dummy2, line
    logical :: exist, parse_is_ok
    real(real64) :: dummyr(7)
 
    push_sub(read_wannier90_files)
 
    w90_num_kpts = product(sys%kpoints%nik_axis(1:3))
    assert(w90_num_kpts == sys%st%nik)
 
 
    w90_num_exclude = 0
 
    ! open nnkp file
    filename = trim(adjustl(w90_prefix)) //'.nnkp'
 
    message(1) = "oct-wannier90: Parsing "//filename
    call messages_info(1)
 
    inquire(file=filename,exist=exist)
    if (.not. exist) then
      message(1) = 'oct-wannier90: Cannot find specified Wannier90 nnkp file.'
      write(message(2),'(a)') 'Please run wannier90.x -pp '// trim(adjustl(w90_prefix)) // ' first.'
      call messages_fatal(2)
    end if
 
    parse_is_ok = .false.
 
    ! check number of k-points
    w90_nnkp = io_open(trim(filename), global_namespace, action='read')
    do
      read(w90_nnkp, *, iostat=io) dummy, dummy1
      if (io == iostat_end) exit
 
      if (dummy == 'begin' .and. dummy1 == 'kpoints') then
        read(w90_nnkp,*) itemp
        if (itemp /= w90_num_kpts) then
          message(1) = 'oct-wannier90: wannier90 setup seems to have been done with a different number of k-points.'
          call messages_fatal(1)
        else
          parse_is_ok = .true.
          exit
        end if
      end if
    end do
    call io_close(w90_nnkp)
 
    if (.not. parse_is_ok) then
      message(1) = 'oct-wannier90: Did not find the kpoints block in nnkp file'
      call messages_fatal(1)
    end if
    parse_is_ok = .false.
 
    ! read from nnkp file
    ! find the nnkpts block
    w90_nnkp = io_open(trim(filename), global_namespace, action='read', position='rewind')
    do
      read(w90_nnkp, *, iostat=io) dummy, dummy1
      if (io  == iostat_end) exit !End of file
 
      if (dummy == 'begin' .and. dummy1 == 'nnkpts') then
        read(w90_nnkp,*) w90_nntot
        safe_allocate(w90_nnk_list(1:5, 1:w90_num_kpts * w90_nntot))
        do ii = 1, w90_num_kpts * w90_nntot
          read(w90_nnkp,*) w90_nnk_list(1:5, ii)
        end do
        !make sure we are at the end of the block
        read(w90_nnkp,*) dummy
        if (dummy /= 'end') then
          message(1) = 'oct-wannier90: There dont seem to be enough k-points in nnkpts file to.'
          call messages_fatal(1)
        end if
        parse_is_ok = .true.
        exit
      end if
    end do
 
    if (.not. parse_is_ok) then
      message(1) = 'oct-wannier90: Did not find nnkpts block in nnkp file'
      call messages_fatal(1)
    end if
 
    ! read from nnkp file
    ! find the exclude block
    safe_allocate(exclude_list(1:sys%st%nst))
    !By default we use all the bands
    exclude_list(1:sys%st%nst) = .false.
    rewind(w90_nnkp)
    do
      read(w90_nnkp, *, iostat=io) dummy, dummy1
      if (io  == iostat_end) exit !End of file
      if (dummy == 'begin' .and. dummy1 == 'exclude_bands') then
        read(w90_nnkp, *) w90_num_exclude
        do ii = 1, w90_num_exclude
          read(w90_nnkp, *) itemp
          if(itemp > sys%st%nst) then
            message(1) = 'oct-wannier90: The exclude_bands list contains a state index higher than the number of states.'
            call messages_fatal(1)
          end if
          exclude_list(itemp) = .true.
        end do
        !make sure we are at the end of the block
        read(w90_nnkp, *) dummy
        if (dummy /= 'end') then
          message(1) = 'oct-wannier90: There dont seem to be enough bands in exclude_bands list.'
          call messages_fatal(1)
        end if
        exit
      end if
    end do
    call io_close(w90_nnkp)
 
    !We get the number of bands
    w90_num_bands = sys%st%nst - w90_num_exclude
 
    safe_allocate(band_index(1:sys%st%nst))
    itemp = 0
    do ii = 1, sys%st%nst
      if (exclude_list(ii)) cycle
      itemp = itemp + 1
      band_index(ii) = itemp
    end do
 
    if (bitand(w90_what, option__wannier90files__w90_amn) /= 0 &
      .or. w90_mode == option__wannier90mode__w90_wannier ) then
      ! parse file again for definitions of projections
      w90_nnkp = io_open(trim(filename), global_namespace, action='read', position='rewind')
 
      do
        read(w90_nnkp, *, iostat=io) dummy, dummy1
        if (io == iostat_end) then !End of file
          message(1) = 'oct-wannier90: Did not find projections block in w90.nnkp file'
          call messages_fatal(1)
        end if
 
        if (dummy == 'begin' .and. (dummy1 == 'projections' .or. dummy1 == 'spinor_projections')) then
 
          if (dummy1 == 'spinor_projections') then
            w90_spinors = .true.
            if (sys%st%d%ispin /= spinors) then
              message(1) = 'oct-wannier90: Spinor = .true. is only valid with spinors wavefunctions.'
              call messages_fatal(1)
            end if
 
            message(1) = 'oct-wannier90: Spinor interface incomplete. Note there is no quantization axis implemented'
            call messages_warning(1)
          else
            if (sys%st%d%ispin == spinors) then
              message(1) = 'oct-wannier90: Octopus has spinors wavefunctions but spinor_projections is not defined.'
              message(2) = 'oct-wannier90: Please check the input file for wannier 90.'
              call messages_fatal(2)
            end if
          end if
 
          read(w90_nnkp, *) w90_nproj
          ! num_wann is given in w90.win, not double checked here
          ! I assume that the wannier90.x -pp run has checked this
          w90_num_wann = w90_nproj
          ! In case of no projections, we use the number of bands
          if(w90_nproj == 0) w90_num_wann = w90_num_bands
 
          safe_allocate(w90_proj_centers(1:3, 1:w90_nproj))
          safe_allocate(w90_proj_lmr(1:w90_nproj, 1:3))
          if (w90_spinors) then
            safe_allocate(w90_spin_proj_component(1:w90_nproj))
          end if
          if (w90_spinors) then
            safe_allocate(w90_spin_proj_axis(1:w90_nproj, 1:3))
          end if
 
          do ii = 1, w90_nproj
            read(w90_nnkp, *) w90_proj_centers(1:3, ii), w90_proj_lmr(ii, 1:3)
            ! skip a line for now
            read(w90_nnkp, *) dummyr(1:7)
            if (w90_spinors) then
              read(w90_nnkp, *) w90_spin_proj_component(ii), w90_spin_proj_axis(ii, 1:3)
              ! use octopus spindim conventions
              if (w90_spin_proj_component(ii) == -1) w90_spin_proj_component(ii) = 2
            end if
          end do
          !make sure we are at the end of the block
          read(w90_nnkp, *) dummy
          if (dummy /= 'end') then
            message(1) = 'oct-wannier90: There dont seem to be enough projections in nnkpts file to.'
            call messages_fatal(1)
          end if
          exit
        end if
      end do
 
      ! look for auto_projection block
      scdm_proj = .false.
      do
        read(w90_nnkp, *, iostat=io) dummy, dummy1
        if (io == iostat_end) exit !End of file
 
        if (dummy == 'begin' .and. dummy1 == 'auto_projections') then
          scdm_proj = .true.
          read(w90_nnkp, *) w90_nproj
          w90_num_wann = w90_nproj
 
          if (.not. w90_scdm) then
            message(1) = 'oct-wannier90: Found auto_projections block. Currently the only implemented automatic way'
            message(2) = 'oct-wannier90: to compute projections is the SCDM method.'
            message(3) = 'oct-wannier90: Please set Wannier90UseSCDM = yes in the inp file.'
            call messages_fatal(3)
          end if
 
          if (w90_nproj /= w90_num_bands) then
            message(1) = 'oct-wannier90: In auto_projections block first row needs to be equal to num_bands.'
            call messages_fatal(1)
          end if
          read(w90_nnkp, *) dummyint
          if (dummyint /= 0) then
            message(1) = 'oct-wannier90: The second row in auto_projections has to be 0, per Wannier90 documentation.'
            call messages_fatal(1)
          end if
        end if
      end do
      call io_close(w90_nnkp)
 
    end if
 
    message(1) = "oct-wannier90: Finished parsing "//filename
    call messages_info(1)
 
    ! Look extra variables variable
    ! open win file
    filename = trim(adjustl(w90_prefix)) //'.win'
    message(1) = "oct-wannier90: Parsing "//filename
    call messages_info(1)
    w90_nnkp = io_open(trim(filename), global_namespace, action='read', position='rewind')
    do
      read(w90_nnkp, fmt='(a)', iostat=io) line
      if (io  == iostat_end) exit !End of file
      if (index(line, '=') > 0) then
        read(line, *, iostat=io) dummy, dummy2, dummy1
      else
        read(line, *, iostat=io) dummy, dummy1
      end if
 
      !Spin
      if (dummy == 'spin') then
        if (sys%st%d%ispin /= spin_polarized) then
          message(1) = 'oct-wannier90: The variable spin is set for a non spin-polarized calculation.'
          call messages_fatal(1)
        end if
 
        if (dummy1 == 'up') then
          spin_channel_win = 1
        else if (dummy1 == 'down') then
          spin_channel_win = 2
        else
          message(1) = 'oct-wannier90: Error parsing the variable spin.'
          call messages_fatal(1)
        end if
        if (spin_channel_win /= w90_spin_channel) then
          message(1) = 'spin polarization input from .win and Wannier90SpinChannel do not agree.'
          call messages_fatal(1)
        end if
      end if
    end do
    call io_close(w90_nnkp)
 
    if (sys%st%d%ispin == spin_polarized) then
      write(message(1), '(a,i1)') 'oct-wannier90: Using spin channel ', w90_spin_channel
      call messages_info(1)
    end if
 
    message(1) = "oct-wannier90: Finished parsing "//filename
    call messages_info(1)
 
    pop_sub(read_wannier90_files)
 
  end subroutine read_wannier90_files
 
  ! --------------------------------------------------------------------------
  subroutine create_wannier90_mmn(mesh, st)
    class(mesh_t),               intent(in) :: mesh
    type(states_elec_t), target, intent(in) :: st
 
    integer ::  ist, jst, ik, ip, w90_mmn, iknn, ib
    real(real64)   ::  Gcart(3)
    integer ::  G(3)
    character(len=80) :: filename
    complex(real64), allocatable :: overlap(:), ss(:)
    complex(real64), allocatable :: psin(:,:), phase(:)
    type(wfs_elec_t), pointer :: batch
    integer :: inode, node_fr, node_to
    type(mpi_request) :: send_req
 
    push_sub(create_wannier90_mmn)
 
    call profiling_in("W90_MMN")
 
    if (st%parallel_in_states) then
      call messages_not_implemented("w90_mmn output with states parallelization")
    end if
 
    message(1) = "Info: Computing the overlap matrix"
    call messages_info(1)
 
 
    filename = './'// trim(adjustl(w90_prefix))//'.mmn'
    w90_mmn = io_open(trim(filename), global_namespace, action='write')
 
    ! write header
    if (st%system_grp%is_root()) then
      write(w90_mmn,*) 'Created by oct-wannier90'
      write(w90_mmn,*) w90_num_bands, w90_num_kpts, w90_nntot
    end if
 
    safe_allocate(psin(1:mesh%np, 1:st%d%dim))
    safe_allocate(phase(1:mesh%np))
    safe_allocate(overlap(1:w90_num_bands))
    safe_allocate(ss(1:st%nst))
 
 
    ! loop over the pairs specified in the nnkp file (read before in init)
    do ii = 1, w90_num_kpts * w90_nntot
      ik = w90_nnk_list(1, ii)
      iknn = w90_nnk_list(2, ii)
      g(1:3) = w90_nnk_list(3:5, ii)
      if (st%system_grp%is_root()) write(w90_mmn, '(I10,2x,I10,2x,I3,2x,I3,2x,I3)') ik, iknn, g
      ! For mixed periodicity we remove the G vectors along the non-periodic directions
      ! as this correspond to the infinite vacuum limit
      g(sys%space%periodic_dim+1:sys%space%dim) = 0
 
      ! For spin-polarized calculations, we select the right k-point
      if (sys%st%d%ispin == spin_polarized) then
        ik = (ik-1)*2 + w90_spin_channel
        iknn = (iknn-1)*2 + w90_spin_channel
      end if
 
 
      ! Only treat local k-points
      if(ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
 
        ! Wannier90 treats everything fully periodic
        ! Conversion is done with the 3D "correct" klattice
        call kpoints_to_absolute(sys%ions%latt, real(g, real64) , gcart)
 
        ! Phase that gives u_{n,k+G}(r) from u_{nk}(r)
        ! Here the minus sign of Octopus cancels with minus sign of the input G
        ! (ik and iknn correspond in the list to minus the k-points in Octopus)
        if (any(g /= 0)) then
          !$omp parallel do
          do ip = 1, mesh%np
            phase(ip) = exp(-m_zi*dot_product(mesh%x(1:3, ip), gcart(1:3)))
          end do
        end if
 
      end if
 
      ! Loop over distributed bands
      ! Slow index in M_{mn}^{k,b}, so \psi_n
      do jst = 1, st%nst
        if (exclude_list(jst)) cycle
 
        ! Communication for the local states
        if ( .not. st%d%kpt%parallel .and. .not. st%parallel_in_states) then
          call states_elec_get_state(st, mesh, jst, iknn, psin)
        else
          node_fr = -1
          node_to = -1
          do inode = 0, st%d%kpt%mpi_grp%size-1
            if(iknn >= st%st_kpt_task(inode,3) .and. iknn <= st%st_kpt_task(inode,4)) then
              node_fr = inode
            end if
            if(ik >= st%st_kpt_task(inode,3) .and. ik <= st%st_kpt_task(inode,4)) then
              node_to = inode
            end if
          end do
          assert(node_fr > -1)
          assert(node_to > -1)
 
          send_req = mpi_request_null
          ! We have locally the k-point
          if (state_kpt_is_local(st, jst, iknn)) then
            call states_elec_get_state(st, mesh, jst, iknn, psin)
            ! We send it only if we don`t want to use it locally
            if(node_to /= st%d%kpt%mpi_grp%rank) then
              call st%d%kpt%mpi_grp%isend(psin, mesh%np*st%d%dim, mpi_double_complex, node_to, send_req)
            end if
          end if
          ! We receive the desired state, only if it is not a local one
          if(node_to == st%d%kpt%mpi_grp%rank .and. node_to /= node_fr) then
            call st%d%kpt%mpi_grp%recv(psin, mesh%np*st%d%dim, mpi_double_complex, node_fr)
          end if
          if (send_req /= mpi_request_null) then
            call st%d%kpt%mpi_grp%wait(send_req)
          end if
        end if
 
        overlap = m_zero
 
        if(ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
 
          ! Do not apply the phase if the phase factor is null
          if (any(g /= 0)) then
            ! add phase
            !$omp parallel
            do idim = 1, st%d%dim
              !$omp do
              do ip = 1, mesh%np
                psin(ip, idim) = psin(ip, idim) * phase(ip)
              end do
            end do
            !$omp end parallel
          end if
 
 
          ! See Eq. (25) in PRB 56, 12847 (1997)
          ! Fast index in M_{mn}^{k,b}, so m.
          do ib = st%group%block_start, st%group%block_end
            batch => st%group%psib(ib, ik)
            if (all(exclude_list(batch%ist(1:batch%nst)))) cycle
            call zmesh_batch_mf_dotp(mesh, batch, psin, ss(1:batch%nst), reduce = .false.)
            do ist = 1, batch%nst
              if (exclude_list(batch%ist(ist))) cycle
              overlap(band_index(batch%ist(ist))) = ss(ist)
            end do
          end do
        end if
 
        call profiling_in("W90_MMN_REDUCE")
        call mesh%allreduce(overlap)
        call profiling_out("W90_MMN_REDUCE")
 
        if(st%d%kpt%parallel) then
          call comm_allreduce(st%d%kpt%mpi_grp, overlap)
        end if
 
        ! write to W90 file
        if (st%system_grp%is_root()) then
          do ist = 1, st%nst
            if (exclude_list(ist)) cycle
            write(w90_mmn,'(e18.10,2x,e18.10)') overlap(band_index(ist))
          end do
        end if
 
      end do !jst
    end do
 
    call io_close(w90_mmn)
 
    safe_deallocate_a(psin)
    safe_deallocate_a(phase)
    safe_deallocate_a(overlap)
    safe_deallocate_a(ss)
 
    call profiling_out("W90_MMN")
 
    pop_sub(create_wannier90_mmn)
 
  end subroutine create_wannier90_mmn
 
  ! --------------------------------------------------------------------------
  subroutine create_wannier90_eig()
    integer ::  ist, ik, w90_eig
    character(len=80) :: filename
 
    push_sub(create_wannier90_eig)
 
    if (sys%st%parallel_in_states) then
      call messages_not_implemented("w90_eig output with states parallelization")
    end if
 
    if (sys%st%system_grp%is_root()) then
      filename = './'//trim(adjustl(w90_prefix))//'.eig'
      w90_eig = io_open(trim(filename), global_namespace, action='write')
      do ik = 1, w90_num_kpts
        do ist = 1, sys%st%nst
          if (exclude_list(ist)) cycle
          if (sys%st%d%ispin /= spin_polarized) then
            write(w90_eig,'(I5,2x,I8,2x,e18.10)') band_index(ist), ik,  &
              units_from_atomic(unit_ev, sys%st%eigenval(ist, ik))
          else
            write(w90_eig,'(I5,2x,I8,2x,e18.10)') band_index(ist), ik,  &
              units_from_atomic(unit_ev, sys%st%eigenval(ist, (ik-1)*2+w90_spin_channel))
          end if
        end do
      end do
 
      call io_close(w90_eig)
    end if
 
    pop_sub(create_wannier90_eig)
  end subroutine create_wannier90_eig
 
  ! --------------------------------------------------------------------------
  subroutine write_unk(space, mesh, st, formatted)
    class(space_t),      intent(in) :: space
    class(mesh_t),       intent(in) :: mesh
    type(states_elec_t), intent(in) :: st
    logical,             intent(in) :: formatted
 
    integer ::  ist, ik, unk_file, ispin
    integer ::  ix, iy, iz
    real(real64) :: w_real, w_imag
    character(len=80) :: filename
    complex(real64), allocatable :: psi(:)
    type(cube_t) :: cube
    type(cube_function_t) :: cf
 
    push_sub(write_unk)
 
    if (st%d%kpt%parallel) then
      call messages_not_implemented("w90_unk output with k-point parallelization")
    end if
 
    if (sys%gr%parallel_in_domains) then
      call messages_not_implemented("w90_unk output with domain parallelization")
    end if
 
    if (st%parallel_in_states) then
      call messages_not_implemented("w90_unk output with states parallelization")
    end if
 
    call messages_experimental("Wannier90Files = w90_unk")
 
 
    safe_allocate(psi(1:mesh%np))
 
    call cube_init(cube, mesh%idx%ll, global_namespace, space, mesh%spacing, &
      mesh%coord_system, need_partition=.not.mesh%parallel_in_domains)
    call cube_init_cube_map(cube, mesh)
 
    call zcube_function_alloc_rs(cube, cf)
 
    do ik = 1, w90_num_kpts
      do ispin = 1, st%d%dim
        if (st%system_grp%is_root()) then
          write(filename, '(a,i5.5,a1,i1)') './UNK', ik,'.', ispin
          ! write header
          if (formatted) then
            unk_file = io_open(trim(filename), global_namespace, action='write', form='formatted')
            write(unk_file, *) mesh%idx%ll(1:mesh%idx%dim), ik, w90_num_bands
          else
            unk_file = io_open(trim(filename), global_namespace, action='write', form='unformatted')
            write(unk_file) mesh%idx%ll(1:mesh%idx%dim), ik, w90_num_bands
          end if
        end if
 
        ! states
        do ist = 1, st%nst
          if (exclude_list(ist)) cycle
 
          if (sys%st%d%ispin /= spin_polarized) then
            call states_elec_get_state(st, mesh, ispin, ist, ik, psi)
          else
            call states_elec_get_state(st, mesh, ispin, ist, (ik-1)*2+w90_spin_channel, psi)
          end if
 
          ! put the density in the cube
          ! Note: At the moment this does not work for domain parallelization
          assert(.not. cube%parallel_in_domains)
          call zmesh_to_cube(mesh, psi, cube, cf)
 
          if (st%system_grp%is_root()) then
            if (formatted) then
              do iz=1,cube%rs_n_global(3)
                do iy = 1,cube%rs_n_global(2)
                  do ix = 1,cube%rs_n_global(1)
                    w_real = real(cf%zrs(ix,iy,iz), real64)
                    w_imag = aimag(cf%zrs(ix,iy,iz))
                    write (unk_file, *) w_real, w_imag
                  end do
                end do
              end do
            else
              write(unk_file) (((cf%zrs(ix,iy,iz), ix=1,cube%rs_n_global(1)), iy=1,cube%rs_n_global(2)), iz=1,cube%rs_n_global(3))
            end if
          end if
        end do
        if (st%system_grp%is_root()) call io_close(unk_file)
      end do
    end do
 
    call zcube_function_free_rs(cube, cf)
    call cube_end(cube)
 
    safe_deallocate_a(psi)
 
    pop_sub(write_unk)
 
  end subroutine write_unk
 
  ! --------------------------------------------------------------------------
  subroutine create_wannier90_amn(space, mesh, latt, st, kpoints)
    class(space_t),          intent(in) :: space
    class(mesh_t),           intent(in) :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    type(states_elec_t),     intent(in) :: st
    type(kpoints_t),         intent(in) :: kpoints
 
    integer ::  ist, ik, w90_amn, idim, iw, ip, ik_real
    real(real64)   ::  center(3),  kpoint(3), threshold
    character(len=80) :: filename
    complex(real64), allocatable :: psi(:,:), phase(:), projection(:)
    real(real64), allocatable ::  ylm(:)
    type(orbitalset_t), allocatable :: orbitals(:)
 
    push_sub(create_wannier90_amn)
    call profiling_in("W90_AMN")
 
    if (st%parallel_in_states) then
      call messages_not_implemented("w90_amn output with states parallelization")
    end if
 
    filename = './'// trim(adjustl(w90_prefix))//'.amn'
    w90_amn = io_open(trim(filename), global_namespace, action='write')
 
    ! write header
    if (st%system_grp%is_root()) then
      write(w90_amn,*) 'Created by oct-wannier90'
      write(w90_amn,*)  w90_num_bands, w90_num_kpts, w90_num_wann
    end if
 
    if (scdm_proj) then
 
      message(1) = "Info: Writing projections obtained from SCDM."
      call messages_info(1)
 
 
      do ik = 1, w90_num_kpts
        do ist = 1, st%nst
          if (exclude_list(ist)) cycle
          if (st%system_grp%is_root()) then
            do iw = 1, w90_nproj
              write (w90_amn,'(I5,2x,I5,2x,I5,2x,e18.10,2x,e18.10)') band_index(ist), iw, ik, uk(band_index(ist),iw,ik)
            end do
          end if
        end do !ist
      end do! ik
 
    else
 
      message(1) = "Info: Computing the projection matrix"
      call messages_info(1)
 
      !We use the variable AOThreshold to deterine the threshold on the radii of the atomic orbitals
      call parse_variable(global_namespace, 'AOThreshold', 0.001_real64, threshold)
 
      safe_allocate(orbitals(1:w90_nproj))
      ! precompute orbitals
      do iw=1, w90_nproj
        call orbitalset_init(orbitals(iw))
 
        orbitals(iw)%norbs = 1
        orbitals(iw)%ndim = 1
        orbitals(iw)%radius = -log(threshold)
        orbitals(iw)%use_submesh = .false.
 
        ! cartesian coordinate of orbital center
        center(1:3) = latt%red_to_cart(w90_proj_centers(1:3, iw))
        call submesh_init(orbitals(iw)%sphere, space, mesh, latt, center, orbitals(iw)%radius)
 
        ! get dorb as submesh points
        safe_allocate(ylm(1:orbitals(iw)%sphere%np))
        ! (this is a routine from pwscf)
        call ylm_wannier(ylm, w90_proj_lmr(iw,1), w90_proj_lmr(iw,2), &
          orbitals(iw)%sphere%r, orbitals(iw)%sphere%rel_x, orbitals(iw)%sphere%np)
 
        ! apply radial function
        if (w90_proj_lmr(iw,3) == 1) then
          do ip = 1,orbitals(iw)%sphere%np
            ylm(ip) = ylm(ip)*m_two*exp(-orbitals(iw)%sphere%r(ip))
          end do
        else
          call messages_not_implemented("oct-wannier90: r/=1 for the radial part")
        end if
 
        safe_allocate(orbitals(iw)%zorb(1:orbitals(iw)%sphere%np, 1, 1))
        orbitals(iw)%zorb(1:orbitals(iw)%sphere%np, 1, 1) = ylm(1:orbitals(iw)%sphere%np)
        safe_deallocate_a(ylm)
 
        safe_allocate(orbitals(iw)%phase(1:orbitals(iw)%sphere%np, st%d%kpt%start:st%d%kpt%end))
        orbitals(iw)%phase(:,:) = m_z0
        safe_allocate(orbitals(iw)%eorb_mesh(1:mesh%np, 1, 1, st%d%kpt%start:st%d%kpt%end))
        orbitals(iw)%eorb_mesh(:,:,:,:) = m_z0
 
        call orbitalset_update_phase(orbitals(iw), space%dim, st%d%kpt, kpoints, st%d%ispin == spin_polarized, &
          kpt_max = w90_num_kpts)
 
      end do
 
      safe_allocate(psi(1:mesh%np, 1:st%d%dim))
      safe_allocate(phase(1:mesh%np))
      safe_allocate(projection(1:w90_nproj))
 
      do ik = 1, w90_num_kpts
        kpoint(1:space%dim) = kpoints%get_point(ik)
        !$omp parallel do
        do ip = 1, mesh%np
          phase(ip) = exp(-m_zi* sum(mesh%x(1:space%dim, ip) * kpoint(1:space%dim)))
        end do
 
        !For spin-polarized calculations, we select the right k-point
        if (st%d%ispin == spin_polarized) then
          ik_real = (ik-1)*2 + w90_spin_channel
        else
          ik_real = ik
        end if
 
 
        do ist = 1, st%nst
          if (exclude_list(ist)) cycle
 
          projection = m_zero
 
          if(ik_real >= st%d%kpt%start .and. ik_real <= st%d%kpt%end) then
            call states_elec_get_state(st, mesh, ist, ik_real, psi)
 
            !$omp parallel
            do idim = 1, st%d%dim
              !The minus sign is here is for the wrong convention of Octopus
              !$omp do
              do ip = 1, mesh%np
                psi(ip, idim) = psi(ip, idim)*phase(ip)
              end do
            end do
            !$omp end parallel
 
            do iw = 1, w90_nproj
              idim = 1
              if (w90_spinors) idim = w90_spin_proj_component(iw)
 
              !At the moment the orbitals do not depend on idim
              !The idim index for eorb_mesh would be for a spin-resolved orbital like j=1/2
              projection(iw) = zmf_dotp(mesh, psi(1:mesh%np,idim), &
                orbitals(iw)%eorb_mesh(1:mesh%np,1,1,ik_real), reduce = .false.)
            end do
 
            call profiling_in("W90_AMN_REDUCE")
            call mesh%allreduce(projection)
            call profiling_out("W90_AMN_REDUCE")
          end if
 
          if(st%d%kpt%parallel) then
            call comm_allreduce(st%d%kpt%mpi_grp, projection)
          end if
 
          if (st%system_grp%is_root()) then
            do iw = 1, w90_nproj
              write (w90_amn,'(I5,2x,I5,2x,I5,2x,e18.10,2x,e18.10)') band_index(ist), iw, ik, projection(iw)
            end do
          end if
        end do !ist
      end do !ik
 
      safe_deallocate_a(psi)
      safe_deallocate_a(phase)
      safe_deallocate_a(projection)
 
      do iw = 1, w90_nproj
        call orbitalset_end(orbitals(iw))
      end do
      safe_deallocate_a(orbitals)
    end if
 
    call io_close(w90_amn)
 
    call profiling_out("W90_AMN")
 
    pop_sub(create_wannier90_amn)
 
  end subroutine create_wannier90_amn
 
  ! --------------------------------------------------------------------------
  subroutine create_wannier90_spn(mesh, st)
    class(mesh_t),               intent(in) :: mesh
    type(states_elec_t), target, intent(in) :: st
 
    integer ::  ist, jst, ik, w90_spn, counter
    character(len=80) :: filename
    complex(real64), allocatable :: spin(:,:,:)
    complex(real64), allocatable :: psim(:,:), psin(:,:)
    complex(real64) :: dot_upup, dot_updown, dot_downup, dot_downdown
 
    push_sub(create_wannier90_spn)
    call profiling_in("W90_SPN")
 
    assert(st%d%ispin == spinors)
 
    if (st%parallel_in_states) then
      call messages_not_implemented("w90_spn output with states parallelization")
    end if
 
    message(1) = "Info: Computing the spin file"
    call messages_info(1)
 
    filename = './'// trim(adjustl(w90_prefix))//'.spn'
    w90_spn = io_open(trim(filename), global_namespace, action='write')
 
    ! write header
    if (st%system_grp%is_root()) then
      write(w90_spn,*) 'Created by oct-wannier90'
      write(w90_spn,*) w90_num_bands, w90_num_kpts
    end if
 
    safe_allocate(psim(1:mesh%np, 1:st%d%dim))
    safe_allocate(psin(1:mesh%np, 1:st%d%dim))
    safe_allocate(spin(1:3, 1:(w90_num_bands*(w90_num_bands+1))/2, 1:w90_num_kpts))
    spin = m_zero
 
    ! loop over the pairs specified in the nnkp file (read before in init)
    do ik = st%d%kpt%start, st%d%kpt%end
      counter = 0
      do jst = 1, st%nst
        if (exclude_list(jst)) cycle
 
        call states_elec_get_state(st, mesh, jst, ik, psim)
        do ist = 1, jst
          if (exclude_list(ist)) cycle
 
          counter = counter + 1
 
          call states_elec_get_state(st, mesh, ist, ik, psin)
 
          dot_upup = zmf_dotp(mesh, psin(:, 1), psim(:, 1), reduce = .false.)
          dot_downdown = zmf_dotp(mesh, psin(:, 2), psim(:, 2), reduce = .false.)
          dot_updown =  zmf_dotp(mesh, psin(:, 1), psim(:, 2), reduce = .false.)
          dot_downup =  zmf_dotp(mesh, psin(:, 2), psim(:, 1), reduce = .false.)
 
          spin(1, counter, ik) = dot_updown + dot_downup
          spin(2, counter, ik) = -m_zi * dot_updown + m_zi * dot_downup
          spin(3, counter, ik) = dot_upup - dot_downdown
        end do !ist
      end do
    end do
 
    call profiling_in("W90_SPN_REDUCE")
    call mesh%allreduce(spin)
 
    if(st%d%kpt%parallel) then
      call comm_allreduce(st%d%kpt%mpi_grp, spin)
    end if
    call profiling_out("W90_SPN_REDUCE")
 
    ! write to W90 file
    if (st%system_grp%is_root()) then
      do ik = 1, w90_num_kpts
        counter = 0
        do jst = 1, st%nst
          if (exclude_list(jst)) cycle
 
          do ist = 1, jst
            if (exclude_list(ist)) cycle
 
            counter = counter + 1
            write(w90_spn, '(e18.10,2x,e18.10)')  spin(1, counter, ik)
            write(w90_spn, '(e18.10,2x,e18.10)')  spin(2, counter, ik)
            write(w90_spn, '(e18.10,2x,e18.10)')  spin(3, counter, ik)
          end do
        end do
      end do !ik
    end if
 
    call io_close(w90_spn)
 
    safe_deallocate_a(psim)
    safe_deallocate_a(psin)
    safe_deallocate_a(spin)
 
    call profiling_out("W90_SPN")
 
    pop_sub(create_wannier90_spn)
  end subroutine create_wannier90_spn
 
 
  ! --------------------------------------------------------------------------
  subroutine generate_wannier_states(space, mesh, ions, st, kpoints)
    class(space_t),         intent(in) :: space
    class(mesh_t),          intent(in) :: mesh
    type(ions_t),           intent(in) :: ions
    type(states_elec_t),    intent(in) :: st
    type(kpoints_t),        intent(in) :: kpoints
 
    integer :: w90_u_mat, w90_xyz, nwann, nik
    integer :: ik, iw, iw2, ip, ipmax, rankmax, idmmax
    real(real64), allocatable :: centers(:,:), supercell_centers(:,:), new_centers(:,:)
    complex(real64), allocatable :: Umnk(:,:,:)
    complex(real64), allocatable :: zwn(:,:,:), psi(:,:), phase(:,:), phase_bloch(:), zwn_bloch(:,:,:)
    character(len=MAX_PATH_LEN) :: fname
    real(real64) :: kpoint(3), wmod, wmodmax, xx(space%dim)
    character(len=2) :: dum
    logical :: exist
    type(unit_t) :: fn_unit
    complex(real64) :: scal
    type(block_t) :: blk
    integer :: supercell(space%dim), ii, jj, kk, ncols, Nreplica, irep, irepmax
    real(real64) :: min_image_displ(3), offset(space%dim)
    integer :: jw
    integer, allocatable :: parent(:)
    real(real64), parameter :: tol_cluster = 0.25_real64 ! In Bohr
 
 
    push_sub(generate_wannier_states)
 
    message(1) = "oct-wannier90: Constructing the Wannier states from the U matrix."
    call messages_info(1)
 
    inquire(file=trim(trim(adjustl(w90_prefix))//'_centres.xyz'),exist=exist)
    if (.not. exist) then
      message(1) = 'oct-wannier90: Cannot find the Wannier90 file seedname_centres.xyz.'
      write(message(2),'(a)') 'Please run wannier90.x with "write_xyz=.true." in '// trim(adjustl(w90_prefix)) // '.'
      call messages_fatal(2)
    end if
 
    w90_xyz = io_open(trim(trim(adjustl(w90_prefix))//'_centres.xyz'), global_namespace, action='read')
 
    safe_allocate(centers(1:3, 1:w90_num_wann))
    !Skip two lines
    read(w90_xyz, *)
    read(w90_xyz, *)
    do iw = 1, w90_num_wann
      read(w90_xyz, *) dum, centers(1:3, iw)
      ! Wannier90 outputs the coordinates in angstrom
      centers(1:3, iw) = units_to_atomic(unit_angstrom, centers(1:3, iw))
    end do
    call io_close(w90_xyz)
 
    ! In order to find "more consistent" Wannier centers, the centers are clustered in groups,
    ! Using periodic boundary condition to find the centers that are the closest
    safe_allocate(parent(1:w90_num_wann))
    do iw = 1, w90_num_wann
      parent(iw) = iw
      do jw = 1, iw-1
        min_image_displ(:) = ions%latt%cart_to_red(centers(:, iw) - centers(:, jw))
        min_image_displ(1:space%periodic_dim) =  min_image_displ(1:space%periodic_dim) - nint(min_image_displ(1:space%periodic_dim))
        min_image_displ = ions%latt%red_to_cart(min_image_displ)
        if (norm2(min_image_displ) < tol_cluster) then
          parent(iw) = parent(jw)
          exit
        end if
      end do
    end do
 
    message(1) = "Info : Clustering of the Wannier centers"
    call messages_info(1)
 
    safe_allocate(new_centers(1:3, 1:w90_num_wann))
    do iw = 1, w90_num_wann
      write(message(1), '(a,i0,a,3(f7.3,a))') 'Info : Original Wannier center ', &
        iw, ' (', centers(1, iw), ', ', centers(2, iw), ', ', centers(3, iw), ')'
 
      if (parent(iw) == iw) then
        ! We are computing |w_{n,0}>, so we need to take only a shift to get to the central cell
        new_centers(1:3, iw) = ions%latt%fold_into_cell(centers(1:3, iw))
      else
        ! This center has a closeby periodic parent. We then move it to be as close as possible to its "parent"
        ! Thanks to the loop above, this should always be inside of very close to a parent inside the unit cell
        min_image_displ(:) = ions%latt%cart_to_red(centers(:, iw) - new_centers(:, parent(iw)))
        min_image_displ(1:space%periodic_dim) =  real(nint(min_image_displ(1:space%periodic_dim)), real64)
        min_image_displ = ions%latt%red_to_cart(min_image_displ)
        new_centers(1:3, iw) = centers(1:3, iw) - min_image_displ
      end if
 
      write(message(2), '(a,i0,a,3(f7.3,a))') 'Info : New Wannier center     ', &
        iw, ' (', new_centers(1, iw), ', ', new_centers(2, iw), ', ', new_centers(3, iw), ')'
      write(message(3), '(a)') ''
      call messages_info(3)
    end do
 
    ! Getting the offset, for the correct phase factor
    centers = new_centers - centers
    safe_deallocate_a(parent)
    safe_deallocate_a(new_centers)
 
    inquire(file=trim(trim(adjustl(w90_prefix))//'_u_dis.mat'),exist=exist)
    if (exist) then
      message(1) = 'oct-wannier90: Calculation of Wannier states with disentanglement is not yet supported.'
      call messages_fatal(1)
    end if
 
    inquire(file=trim(trim(adjustl(w90_prefix))//'_u.mat'),exist=exist)
    if (.not. exist) then
      message(1) = 'oct-wannier90: Cannot find the Wannier90 seedname_u.mat file.'
      write(message(2),'(a)') 'Please run wannier90.x with "write_u_matrices=.true." in '// trim(adjustl(w90_prefix)) // '.'
      call messages_fatal(2)
    end if
    w90_u_mat = io_open(trim(trim(adjustl(w90_prefix))//'_u.mat'), global_namespace, action='read')
 
 
    !Skip one line
    read(w90_u_mat, *)
    !Read num_kpts, num_wann, num_wann for consistency check
    read(w90_u_mat, *) nik, nwann, nwann
    if (nik /= w90_num_kpts .or. nwann /= w90_num_wann) then
      message(1) = "The file contains U matrices is inconsistent with the .win file."
      call messages_fatal(1)
    end if
 
    safe_allocate(umnk(1:w90_num_wann, 1:w90_num_wann, 1:w90_num_kpts))
 
    do ik = 1, w90_num_kpts
      !Skip one line
      read(w90_u_mat, *)
      !Skip one line (k-point coordinate)
      read(w90_u_mat, *)
      read(w90_u_mat, '(f15.10,sp,f15.10)') ((umnk(iw, iw2, ik), iw=1, w90_num_wann), iw2=1, w90_num_wann)
    end do
 
    call io_close(w90_u_mat)
 
    !We read the output format for the Wannier states
    call parse_variable(global_namespace, 'OutputFormat', 0, how)
    if (how == 0) then
      message(1) = "OutputFormat must be specified for outputing Wannier functions."
      call messages_fatal(1)
    end if
 
    !%Variable Wannier90Supercell
    !%Type block
    !%Section Utilities::oct-wannier90
    !%Description
    !% This block allows to specify the size of the supercell used to compute the Wannier functions
    !% If not specified, the code uses a default 1x1x1 cell, i.e., it plots the Wannier90 in the primitive cell.
    !%End
    if (parse_is_defined(sys%namespace, 'Wannier90Supercell')) then
      if (parse_block(sys%namespace, 'Wannier90Supercell', blk) == 0) then
        ncols = parse_block_cols(blk, 0)
        if (ncols /= space%dim) then
          write(message(1),'(a,i3,a,i3)') 'Wannier90Supercell has ', ncols, ' columns but must have ', sys%space%dim
          call messages_fatal(1, namespace=sys%namespace)
        end if
        do ii = 1, space%dim
          call parse_block_integer(blk, 0, ii - 1, supercell(ii))
        end do
 
        call parse_block_end(blk)
      end if
    else
      supercell = 1
    end if
 
    nreplica = product(supercell)
 
    ! The center of each replica of the unit cell
    safe_allocate(supercell_centers(1:space%dim, 1:nreplica))
    offset(1:space%dim) = -floor((real(supercell(1:space%dim), real64) - m_one) / m_two)
    irep = 1
    do ii = 0, supercell(1)-1
      do jj = 0, supercell(2)-1
        do kk = 0, supercell(3)-1
          supercell_centers(:, irep) = ions%latt%red_to_cart(offset &
            + [real(ii, real64), real(jj, real64), real(kk, real64)])
          irep = irep + 1
        end do
      end do
    end do
 
 
 
    call io_mkdir('wannier', global_namespace)
 
    !Computing the Wannier states in the primitive cell, from the U matrices
    safe_allocate(zwn(1:mesh%np, 1:nreplica, 1:st%d%dim))
    safe_allocate(psi(1:mesh%np, 1:st%d%dim))
    safe_allocate(phase(1:mesh%np, 1:nreplica))
    safe_allocate(phase_bloch(1:mesh%np))
    if (w90_bloch_sums) then
      safe_allocate(zwn_bloch(1:mesh%np, 1:st%d%dim, st%d%kpt%start:st%d%kpt%end))
    end if
 
    do iw = 1, w90_num_wann
 
      zwn(:, :, :) = m_z0
 
      do ik = 1, w90_num_kpts
 
        if (.not. (ik >= st%d%kpt%start .and. ik <= st%d%kpt%end)) cycle
 
        if (w90_bloch_sums) zwn_bloch(:,:,ik) = m_z0
 
        kpoint(1:space%dim) = kpoints%get_point(ik, absolute_coordinates=.true.)
 
        ! We compute the Wannier orbital on a grid centered around the Wannier function
        ! The minus sign is here is for the wrong convention of Octopus
        do irep = 1, nreplica
          !$omp parallel do
          do ip = 1, mesh%np
            xx = mesh%x(1:space%dim, ip)-centers(1:space%dim, iw) + supercell_centers(1:space%dim, irep)
            phase(ip, irep) = exp(-m_zi* sum( xx * kpoint(1:space%dim)))
          end do
        end do
        if (w90_bloch_sums) then
          do ip = 1, mesh%np
            xx = mesh%x(1:space%dim, ip)
            phase_bloch(ip) = exp(-m_zi* sum( xx * kpoint(1:space%dim)))
          end do
        end if
 
        do iw2 = 1, st%nst
          if (exclude_list(iw2)) cycle
 
          if (st%d%ispin /= spin_polarized) then
            call states_elec_get_state(st, mesh, iw2, ik, psi)
          else
            call states_elec_get_state(st, mesh, iw2, (ik-1)*2+w90_spin_channel, psi)
          end if
 
          !$omp parallel
          do idim = 1, st%d%dim
            do irep = 1, nreplica
              !$omp do
              do ip = 1, mesh%np
                zwn(ip, irep, idim) = zwn(ip, irep, idim) + umnk(band_index(iw2), iw, ik) * psi(ip, idim) * phase(ip, irep)
              end do
            end do
          end do
          !$omp end parallel
 
          ! See Eq. 19 in arXiv:0605539
          if (w90_bloch_sums) then
            !$omp parallel
            do idim = 1, st%d%dim
              !$omp do
              do ip = 1, mesh%np
                zwn_bloch(ip, idim, ik) = zwn_bloch(ip, idim, ik) &
                  + umnk(band_index(iw2), iw, ik) * psi(ip, idim) * phase_bloch(ip)
              end do
            end do
            !$omp end parallel
          end if
        end do!iw2
      end do!ik
 
      if(st%d%kpt%parallel) then
        call comm_allreduce(st%d%kpt%mpi_grp, zwn)
      end if
 
      ! Following what Wannier90 is doing, we fix the global phase by setting the max to be real
      ! We also normalize to the number of k-point at this step
      if (sys%st%d%ispin /= spinors) then
        ipmax = 0
        irepmax = 0
        wmodmax = m_zero
        do irep = 1, nreplica
          do ip = 1, mesh%np
            wmod = real(zwn(ip, irep, 1)*conjg(zwn(ip, irep, 1)), real64)
            if (wmod > wmodmax) then
              ipmax = ip
              irepmax = irep
              wmodmax = wmod
            end if
          end do
        end do
        scal = sqrt(wmodmax)/zwn(ipmax, irepmax, 1)/w90_num_kpts
        call mesh_minmaxloc(mesh, wmodmax, rankmax, mpi_maxloc)
        call mesh%mpi_grp%bcast(scal, 1, mpi_double_complex, rankmax)
        call lalg_scal(mesh%np, nreplica, scal, zwn(:,:,1))
 
        if (w90_bloch_sums) then
          do ik = st%d%kpt%start, st%d%kpt%end
            ipmax = 0
            idmmax= 0
            wmodmax = m_zero
            do idim = 1, st%d%dim
              do ip = 1, mesh%np
                wmod = real(zwn_bloch(ip, idim, ik)*conjg(zwn_bloch(ip, idim, ik)), real64)
                if (wmod > wmodmax) then
                  ipmax = ip
                  wmodmax = wmod
                  idmmax = idim
                end if
              end do
            end do
            scal = sqrt(wmodmax)/zwn_bloch(ipmax, idmmax, ik)/w90_num_kpts
            call mesh_minmaxloc(mesh, wmodmax, rankmax, mpi_maxloc)
            call mesh%mpi_grp%bcast(scal, 1, mpi_double_complex, rankmax)
            call lalg_scal(mesh%np, st%d%dim, scal, zwn_bloch(:,:,ik))
          end do
        end if
      end if
 
      ! Output the Wannier states
      fn_unit = sqrt(units_out%length**(-space%dim))
      do idim = 1, st%d%dim
        if (st%d%ispin == spinors) then
          write(fname, '(a,i3.3,a4,i1)') 'wannier-', iw, '-isp', idim
        else
          write(fname, '(a,i3.3,a4,i1)') 'wannier-', iw
        end if
        if (any(supercell_centers > 1)) then
          call io_function_output_supercell(how, "wannier", trim(fname), mesh, &
            space, ions%latt, zwn(:, :, idim), supercell_centers, supercell, fn_unit, ierr, global_namespace,       &
            pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
 
        else
          call zio_function_output(how, "wannier", trim(fname), global_namespace, space, mesh, &
            zwn(:, 1, idim), fn_unit, ierr, pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
        end if
      end do
 
      ! Output Bloch sums of Wannier states
      if (w90_bloch_sums) then
        do ik = st%d%kpt%start, st%d%kpt%end
          do idim = 1, st%d%dim
            if (st%d%ispin == spinors) then
              write(fname, '(a,i3.3,a4,i1,a3,i5.5)') 'wannier_bloch-', iw, '-isp', idim, '-ik', ik
            else
              write(fname, '(a,i3.3,a4,i1,a3,i5.5)') 'wannier_bloch-', iw, '-ik', ik
            end if
            call zio_function_output(how, "wannier", trim(fname), global_namespace, space, mesh, &
              zwn_bloch(:, idim, ik), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
          end do
        end do
      end if
 
      ! Checking the ratio imag/real
      if (sys%st%d%ispin /= spinors) then
        wmodmax = m_zero
        do irep = 1, nreplica
          do ip = 1, mesh%np
            if(abs(real(zwn(ip, irep, 1), real64)) >= 1e-2_real64) then
              wmodmax = max(wmodmax, abs(aimag(zwn(ip, irep, 1)))/abs(real(zwn(ip, irep, 1), real64)))
            end if
          end do
        end do
        call mesh_minmaxloc(mesh, wmodmax, rankmax, mpi_maxloc)
 
        write(message(1), '(a,i4,a,f11.6)') 'oct-wannier90: Wannier function ', iw, ' Max. Im/Re Ratio = ', wmodmax
        call messages_info(1)
      else
        write(message(1), '(a,i4)') 'oct-wannier90: Wannier function done ', iw
        call messages_info(1)
      end if
    end do
 
    safe_deallocate_a(umnk)
    safe_deallocate_a(zwn)
    safe_deallocate_a(psi)
    safe_deallocate_a(phase)
    safe_deallocate_a(phase_bloch)
    safe_deallocate_a(centers)
    safe_deallocate_a(supercell_centers)
 
    pop_sub(generate_wannier_states)
  end subroutine generate_wannier_states
 
end program wannier90_interface
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: