wannier90lib_interface.F90

Go to the documentation of this file.

!! Copyright (C) 2026 A. Buccheri, J. Reimann
!!
!! 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"
 
 
module wannier90lib_interface_oct_m
  use debug_oct_m
  use electron_space_oct_m,  only: spin_polarized
  use global_oct_m
  use ions_oct_m,            only: ions_t
  use io_oct_m,              only: io_open, io_close
  use kpoints_oct_m,         only: kpoints_t
  use lattice_vectors_oct_m, only: lattice_vectors_t
  use mesh_oct_m,            only: mesh_t
  use messages_oct_m
  use namespace_oct_m,       only: namespace_t
  use profiling_oct_m
  use space_oct_m,           only: space_t
  use states_elec_oct_m,     only: states_elec_t
  use unit_oct_m,            only: units_from_atomic, units_to_atomic
  use unit_system_oct_m,     only: unit_angstrom, unit_ev
  use wannier_calc_oct_m, &
    only: &
    wannier_calc_proj_t, &
    wannier_calc_create_amn, &
    wannier_calc_create_mmn
  use wannier_opts_oct_m,   only: wannier_opts_t
 
  ! External library
  ! Wannier90 lib type and routines are not dressed with the preprocessor guard as they
  ! are limited to this scope, and this module should only get included in compilation
  ! if Wannier90 lib is linked to
  use w90_library
 
  implicit none
  private
  public :: wannier90lib_init_w90main,  &
    wannier90lib_initialize_matrix_data, &
    wannier90lib_create_wannier90_mmn,  &
    wannier90lib_create_wannier90_amn,  &
    wannier90lib_wannierization
 
contains
 
  subroutine wannier90lib_init_w90main(namespace, ions, latt, kpoints, st, &
    inp_options, w90main, stdout, stderr, ierr)
    type(namespace_t),       intent(in) :: namespace
    type(ions_t),            intent(in) :: ions
    type(lattice_vectors_t), intent(in) :: latt  ! will differ from ions%latt by 2*L for non-periodic systems
    type(kpoints_t),         intent(in) :: kpoints
    type(states_elec_t),     intent(in) :: st
    type(wannier_opts_t),    intent(in) :: inp_options
    type(lib_common_type),   intent(inout) :: w90main
    integer,                 intent(inout) :: stdout, stderr, ierr
 
    integer, allocatable :: distk(:)
    integer :: ik, iq, nk
 
    push_sub(wannier90lib_init_w90main)
 
    ! Set required options
    ! Note: num_bands does not include excluded bands!
    ! The order is nst > num_bands > num_wann
    ! TODO: Correctly handle this depending on inputs
    ! TODO: Distinguish everywhere between nst and num_bands
    call w90_set_option(w90main, 'kpoints', -kpoints%reduced%red_point(1:3,:)) ! minus sign for wrong octopus convention
    call w90_set_option(w90main, 'mp_grid', kpoints%nik_axis(1:3))
    call w90_set_option(w90main, 'num_bands', inp_options%num_bands)
    call w90_set_option(w90main, 'num_wann', inp_options%num_wann)
    call w90_set_option(w90main, 'unit_cell_cart', units_from_atomic(unit_angstrom, latt%rlattice(1:3,:)))
 
    nk = kpoints%full%npoints
    safe_allocate(distk(1:nk))
    if (st%d%kpt%parallel) then
      ! pass octopus communicator to wannier90 library
      call w90_set_comm(w90main, st%d%kpt%mpi_grp%comm%MPI_VAL)
      ! Feed kpoint distribution to wannier90 library
      distk = -1
      ! (k, spin)
      do iq = 1, st%nik
        ! k-only
        ik = st%d%get_kpoint_index(iq)
        if (distk(ik) == -1) then
          distk(ik) = st%d%kpt%process(iq)
        else if (distk(ik) /= st%d%kpt%process(iq)) then
          message(1) = 'Error: Wannier90 k-point parallelization: same k-point' // &
            'assigned to different nodes.'
          call messages_fatal(1, namespace=namespace)
        end if
      end do
    else
      ! octopus default communicator is not acceptable in that case
      call w90_set_comm(w90main, st%system_grp%comm%MPI_VAL)
      ! No parallelization, assign all k-points to node 0
      distk = 0
    endif
    call w90_set_option(w90main, 'distk', distk)
    safe_deallocate_a(distk)
 
    ! Set atom data for projection parsing
    call wannier90lib_set_atom_data(w90main, ions, latt)
 
    ! Parse projections from input file
    call wannier90lib_set_projections(w90main, inp_options)
 
    ! Option to dump files
    if (inp_options%dump_inputs) then
      call w90_set_option(w90main, "dump_inputs", .true.)
    end if
 
    ! Initialize library type with seedname
    call w90_input_setopt(w90main, inp_options%prefix, stdout, stderr, ierr)
 
    ! Read optional options from file
    ! Note: Mandatory options are only read from octopus input.
    ! If any mandatory options are set in `.win`, they will be ignored triggering a warning.
    call wannier90lib_warn_inputs(inp_options)
    call w90_input_reader(w90main, stdout, stderr, ierr)
    if (ierr /= 0) then
      write(message(1), '(a,i0)') 'Error in w90_input_reader: ', ierr
      call messages_fatal(1)
    end if
 
    pop_sub(wannier90lib_init_w90main)
  end subroutine wannier90lib_init_w90main
 
 
  subroutine wannier90lib_set_atom_data(w90main, ions, latt)
    type(lib_common_type), intent(inout) :: w90main
    type(ions_t),          intent(in)    :: ions
    type(lattice_vectors_t), intent(in)  :: latt    ! will differ from ions%latt by 2*L for non-periodic systems
 
    integer :: ia, n_atoms
    real(real64), allocatable :: atom_pos_frac(:,:)
    character(len=2), allocatable :: atom_labels(:)
 
    push_sub(wannier90lib_set_atom_data)
 
    n_atoms = ions%natoms
 
    ! Allocate arrays for atom data
    safe_allocate(atom_pos_frac(1:3, 1:n_atoms))
    safe_allocate(atom_labels(1:n_atoms))
 
    ! Extract atomic positions and labels from ions
    do ia = 1, n_atoms
      ! Get Cartesian positions from octopus
      ! Convert to fractional (reduced) coordinates for wannier90
      atom_pos_frac(1:3, ia) = latt%cart_to_red(ions%pos(1:3, ia))
      ! Get atomic labels
      atom_labels(ia) = trim(ions%atom(ia)%label)
    end do
 
    ! Set atom positions in fractional coordinates
    call w90_set_option(w90main, 'atoms_frac', atom_pos_frac)
    ! Set symbols
    call w90_set_option(w90main, 'symbols', atom_labels)
 
    safe_deallocate_a(atom_pos_frac)
    safe_deallocate_a(atom_labels)
 
    pop_sub(wannier90lib_set_atom_data)
  end subroutine wannier90lib_set_atom_data
 
 
  subroutine wannier90lib_set_projections(w90main, inp_options)
    type(lib_common_type), intent(inout) :: w90main
    type(wannier_opts_t),     intent(in)    :: inp_options
 
    character(len=80)       :: filename, line
    integer                 :: w90_win, io
    logical                 :: exist
 
    push_sub(wannier90lib_set_projections)
 
    ! Correctly handle projections by parsing them from .win file
    filename = trim(adjustl(inp_options%prefix)) //'.win'
 
    message(1) = "oct-wannier90: Parsing projections from "//filename
    call messages_info(1)
 
    inquire(file=filename,exist=exist)
    if (.not. exist) then
      message(1) = 'oct-wannier90: Cannot find specified Wannier90 win file.'
      call messages_fatal(1)
    end if
 
    ! read projections and set them in w90main
    w90_win = io_open(trim(filename), action='read')
    do
      read(w90_win, '(A)', iostat=io) line
      if (is_iostat_end(io)) exit
 
      if (trim(line) == 'begin projections') then
        do
          read(w90_win, '(A)') line
          if (trim(line) == 'end projections') exit
          call w90_set_option(w90main, 'projections', trim(line))
        end do
      end if
    end do
    call io_close(w90_win)
 
    pop_sub(wannier90lib_set_projections)
  end subroutine wannier90lib_set_projections
 
  subroutine wannier90lib_initialize_matrix_data(w90main, u_matrix, &
    u_dis_matrix, m_matrix, centers)
    type(lib_common_type),       intent(inout) :: w90main
    complex(real64), contiguous, intent(inout) :: u_matrix(:, :, :)
    complex(real64), contiguous, intent(inout) :: u_dis_matrix(:, :, :)
    complex(real64), contiguous, intent(inout) :: m_matrix(:, :, :, :)
    real(real64), contiguous,    intent(in)    :: centers(:, :)
 
    push_sub(wannier90lib_initialize_matrix_data)
 
    call w90_set_u_matrix(w90main, u_matrix)
    call w90_set_u_opt(w90main, u_dis_matrix)
    call w90_set_m_local(w90main, m_matrix)
    w90main%wannier_data%centres = centers
 
    pop_sub(wannier90lib_initialize_matrix_data)
 
  end subroutine wannier90lib_initialize_matrix_data
 
  subroutine wannier90lib_warn_inputs(inp_options)
    type(wannier_opts_t),     intent(in)    :: inp_options
 
    character(len=80)       :: filename, line, dummy, dummy1
    integer                 :: w90_win, io
    logical                 :: exist
 
    push_sub(wannier90lib_warn_inputs)
 
    filename = trim(adjustl(inp_options%prefix)) //'.win'
 
    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 win file.'
      call messages_fatal(1)
    end if
 
    ! warn which options are not considered from .win file
    w90_win = io_open(trim(filename), action='read')
    do
      read(w90_win, '(A)', iostat=io) line
      if (is_iostat_end(io)) exit
 
      if (trim(line) == 'begin kpoints') then
        message(1) = 'oct-wannier90: Warning: kpoints specified in .win file' // &
          ' are ignored. Using kpoints from octopus input instead.'
        call messages_warning(1)
      else if (trim(line) == 'begin atoms_frac') then
        message(1) = 'oct-wannier90: Warning: atomic positions specified in .win file' // &
          ' are ignored. Using positions from octopus input instead.'
        call messages_warning(1)
      else if (trim(line) == 'begin unit_cell_cart') then
        message(1) = 'oct-wannier90: Warning: unit cell vectors specified in .win file' // &
          ' are ignored. Using unit cell vectors from octopus input instead.'
        call messages_warning(1)
      end if
      if (index(line, '=') <= 0) then
        read(line, *, iostat=io) dummy, dummy1
        if(dummy == 'mp_grid') then
          message(1) = 'oct-wannier90: Warning: mp_grid specified in .win file' // &
            ' is ignored. Using k-point mesh from octopus input instead.'
          call messages_warning(1)
        end if
      end if
    end do
    call io_close(w90_win)
 
    pop_sub(wannier90lib_warn_inputs)
  end subroutine wannier90lib_warn_inputs
 
  subroutine wannier90lib_create_wannier90_amn(w90main, inp_options, &
    band_index, space, mesh, latt, st, kpoints, projection)
    type(lib_common_type),   intent(in)    :: w90main
    type(wannier_opts_t),    intent(in)    :: inp_options
    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, contiguous,     intent(in)    :: band_index(:)
 
    complex(real64), contiguous, intent(out) :: projection(:, :, :)
    type(wannier_calc_proj_t), allocatable   :: proj_input_oct(:)
 
    integer ::  iw
 
    push_sub(wannier90lib_create_wannier90_amn)
    call profiling_in("W90_AMN")
 
    ! Note: in w90main, proj_input contains the projections from input file,
    ! proj contains the projections that are actually used in wannier90
    ! We use proj_input since that is also what is written to file in -pp mode
 
    
    safe_allocate(proj_input_oct(1:w90main%num_proj))
    do iw = 1, w90main%num_proj
      proj_input_oct(iw)%l = w90main%proj_input(iw)%l
      proj_input_oct(iw)%m = w90main%proj_input(iw)%m
      proj_input_oct(iw)%s = w90main%proj_input(iw)%s
      proj_input_oct(iw)%radial = w90main%proj_input(iw)%radial
      proj_input_oct(iw)%site = w90main%proj_input(iw)%site
      proj_input_oct(iw)%s_qaxis = w90main%proj_input(iw)%s_qaxis
      proj_input_oct(iw)%z = w90main%proj_input(iw)%z
      proj_input_oct(iw)%x = w90main%proj_input(iw)%x
      proj_input_oct(iw)%zona = w90main%proj_input(iw)%zona
    end do
 
    call wannier_calc_create_amn(inp_options, &
      band_index, space, mesh, latt, st, kpoints, &
      w90main%wvfn_read%spin_channel, w90main%num_proj, &
      proj_input_oct, projection)
 
    safe_deallocate_a(proj_input_oct)
 
    call profiling_out("W90_AMN")
 
    pop_sub(wannier90lib_create_wannier90_amn)
 
  end subroutine wannier90lib_create_wannier90_amn
 
  subroutine wannier90lib_create_wannier90_mmn(w90main, band_index, mesh, st, latt, overlap)
    type(lib_common_type),       intent(in) :: w90main
    integer, contiguous,         intent(in) :: band_index(:)
    class(mesh_t),               intent(in) :: mesh
    type(states_elec_t), target, intent(in) :: st
    type(lattice_vectors_t),     intent(in) :: latt    ! will differ from ions%latt by 2*L for non-periodic systems
 
    complex(real64), contiguous, intent(out) :: overlap(:, :, :, :)
 
    push_sub(wannier90lib_create_wannier90_mmn)
 
    call wannier_calc_create_mmn(band_index, &
      mesh, st, latt, w90main%wvfn_read%spin_channel, &
      w90main%kmesh_info%nnlist, w90main%kmesh_info%nncell, &
      overlap)
 
    pop_sub(wannier90lib_create_wannier90_mmn)
 
  end subroutine wannier90lib_create_wannier90_mmn
 
  subroutine wannier90lib_wannierization(w90main, gr, space, st, latt, kpoints, options, band_index, &
    m_matrix, u_matrix_opt, last_u_matrix, warm_start)
    type(lib_common_type),       intent(inout) :: w90main
    class(mesh_t),               intent(in)    :: gr
    class(space_t),              intent(in)    :: space
    type(states_elec_t),         intent(in)    :: st
    type(lattice_vectors_t),     intent(in)    :: latt  ! will differ from ions%latt by 2*L for non-periodic systems
    type(kpoints_t),             intent(in)    :: kpoints
    type(wannier_opts_t),        intent(in)    :: options
    integer, contiguous,         intent(in)    :: band_index(:)
    complex(real64), contiguous, intent(inout) :: m_matrix(:, :, :, :)
    complex(real64), contiguous, intent(inout) :: u_matrix_opt(:, :, :)
    complex(real64), contiguous, intent(in)    :: last_u_matrix(:, :, :)
    logical,          optional,  intent(in)    :: warm_start
 
    integer                     :: ierr
    integer                     :: null_stdout
    logical                     :: warm_start_
 
    push_sub(wannier90lib_wannierization)
 
    warm_start_ = .true.
    if (present(warm_start)) warm_start_ = warm_start
 
    ! Calls to libwannier90 produce output on every call. We suppress this output here
    null_stdout = io_open('/dev/null', action='write', die=.false.)
 
    ! Compute quantities required by wannier90
    ! m matrix is always needed, as the overlaps change with time
    call wannier90lib_create_wannier90_mmn(w90main, band_index, gr, st, latt, m_matrix)
 
    ! Perform wannierization
    if (.not. warm_start_) then
      ! Cold start: initialize U from AMN projections
      call wannier90lib_create_wannier90_amn(w90main, options, band_index, space, &
        gr, latt, st, kpoints, u_matrix_opt)
    else
      ! Warm start: initialize U from previous wannierization result
      u_matrix_opt(:, :, :) = last_u_matrix(:, :, :)
    end if
    ! Not needed for now, since num_bands \= num_wann is not yet supported
    ! call w90_disentangle(w90main, null_stdout, stderr, ierr)
    ! if (ierr /= 0) then
    !   write(message(1), '(a,i0)') 'Error in w90_disentangle: ', ierr
    !   call messages_fatal(1)
    ! end if
    call w90_project_overlap(w90main, null_stdout, stderr, ierr) 
    if (ierr /= 0) then
      write(message(1), '(a,i0)') 'Error in w90_project_overlap: ', ierr
      call messages_fatal(1)
    end if
    call w90_wannierise(w90main, null_stdout, stderr, ierr)
    if (ierr /= 0) then
      write(message(1), '(a,i0)') 'Error in w90_wannierise: ', ierr
      call messages_fatal(1)
    end if
 
    call io_close(null_stdout)
 
    pop_sub(wannier90lib_wannierization)
  end subroutine wannier90lib_wannierization
 
end module wannier90lib_interface_oct_m