plot_wannier Subroutine

  subroutine plot_wannier(wannier_plot, wvfn_read, wannier_data, print_output, u_matrix_opt, &
                          dis_manifold, real_lattice, atom_data, kpt_latt, u_matrix, num_kpts, &
                          num_bands, num_wann, have_disentangled, spinors, bohr, stdout, seedname, &
                          timer, dist_k, error, comm)
    !================================================!
    !! Plot the WF in Xcrysden format
    !! based on code written by Michel Posternak
    !
    !================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_i, twopi, cmplx_1
    use w90_io, only: io_date, io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: wannier_data_type, atom_data_type, dis_manifold_type, print_output_type, &
                         timer_list_type
    use w90_wannier90_types, only: wvfn_read_type, wannier_plot_type
    use w90_comms, only: w90_comm_type
    use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_file, &
                         set_error_warn

    implicit none

    ! arguments
    type(atom_data_type), intent(in) :: atom_data
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(print_output_type), intent(in) :: print_output
    type(wannier_data_type), intent(in) :: wannier_data
    type(wannier_plot_type), intent(in) :: wannier_plot
    type(wvfn_read_type), intent(in) :: wvfn_read
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    complex(kind=dp), intent(in) :: u_matrix(:, :, :)
    complex(kind=dp), intent(in) :: u_matrix_opt(:, :, :)

    real(kind=dp), intent(in) :: bohr
    real(kind=dp), intent(in) :: kpt_latt(:, :)
    real(kind=dp), intent(in) :: real_lattice(3, 3)

    integer, intent(in) :: dist_k(:)
    integer, intent(in) :: num_bands
    integer, intent(in) :: num_kpts
    integer, intent(in) :: num_wann
    integer, intent(in) :: stdout

    logical, intent(in) :: have_disentangled
    logical, intent(in) :: spinors

    character(len=50), intent(in)  :: seedname

    ! local variables
    real(kind=dp) :: tmax, tmaxx, x_0ang, y_0ang, z_0ang
    real(kind=dp) :: fxcry(3), dirl(3, 3), ratmax, ratio
    real(kind=dp) :: upspinor, dnspinor, upphase, dnphase

    complex(kind=dp), allocatable :: wann_func(:, :, :, :)
    complex(kind=dp), allocatable :: r_wvfn(:, :)
    complex(kind=dp), allocatable :: r_wvfn_tmp(:, :)
    complex(kind=dp), allocatable :: wann_func_nc(:, :, :, :, :) ! add the spinor dim.
    complex(kind=dp), allocatable :: r_wvfn_nc(:, :, :) ! add the spinor dim.
    complex(kind=dp), allocatable :: r_wvfn_tmp_nc(:, :, :) ! add the spinor dim.
    complex(kind=dp) :: catmp, wmod

    logical :: have_file, on_root

    integer :: num_nodes, my_node_id
    integer :: i, j, nsp, nat, nbnd, ierr
    integer :: loop_kpt, nk, ngx, ngy, ngz, nxx, nyy, nzz
    integer :: loop_b, nx, ny, nz, npoint, file_unit, loop_w, num_inc
    integer :: wann_plot_num

    character(len=11) :: wfnname
    character(len=60) :: wanxsf, wancube
    character(len=9)  :: cdate, ctime
    logical           :: inc_band(num_bands)

    ! Pre-contracted wavefunction: band sum done once per unit cell point
    complex(kind=dp), allocatable :: c_wvfn(:, :)
    complex(kind=dp), allocatable :: c_wvfn_nc(:, :, :)
    complex(kind=dp) :: uw, phase_yz
    integer :: ngrid

    ! Factored phase arrays
    complex(kind=dp), allocatable :: phase_x(:), phase_y(:), phase_z(:)
    integer :: nxx_lo, nxx_hi, nyy_lo, nyy_hi, nzz_lo, nzz_hi

    num_nodes = mpisize(comm)
    my_node_id = mpirank(comm)

    on_root = .false.
    if (my_node_id == 0) on_root = .true.

    !
    if (print_output%timing_level > 1) call io_stopwatch_start('plot: wannier', timer)
    !
    associate (ngs => wannier_plot%supercell)
      !
      if (.not. spinors) then
        write (wfnname, 200) 1, wvfn_read%spin_channel
      else
        write (wfnname, 199) 1
      end if
      inquire (file=wfnname, exist=have_file)
      if (.not. have_file) then
        call set_error_file(error, 'plot_wannier: file '//wfnname//' not found', comm)
        return
      end if

      if (wvfn_read%formatted) then
        open (newunit=file_unit, file=wfnname, form='formatted', status='old', iostat=ierr)
      else
        open (newunit=file_unit, file=wfnname, form='unformatted', status='old', iostat=ierr)
      end if
      if (ierr /= 0) then
        call set_error_file(error, 'plot_wannier: could not open file '//wfnname, comm)
        return
      end if
      if (wvfn_read%formatted) then
        read (file_unit, *, iostat=ierr) ngx, ngy, ngz, nk, nbnd
      else
        read (file_unit, iostat=ierr) ngx, ngy, ngz, nk, nbnd
      end if
      close (file_unit)
      if (ierr /= 0) then
        call set_error_file(error, 'plot_wannier: error reading file '//wfnname, comm)
        return
      end if

200   format('UNK', i5.5, '.', i1)
199   format('UNK', i5.5, '.', 'NC')

      if (allocated(wannier_plot%list)) then
        wann_plot_num = size(wannier_plot%list)
      else
        wann_plot_num = 0
      end if

      ! Supercell grid bounds
      nxx_lo = -((ngs(1))/2)*ngx
      nxx_hi = ((ngs(1) + 1)/2)*ngx - 1
      nyy_lo = -((ngs(2))/2)*ngy
      nyy_hi = ((ngs(2) + 1)/2)*ngy - 1
      nzz_lo = -((ngs(3))/2)*ngz
      nzz_hi = ((ngs(3) + 1)/2)*ngz - 1
      ngrid = ngx*ngy*ngz

      allocate (wann_func(nxx_lo:nxx_hi, nyy_lo:nyy_hi, nzz_lo:nzz_hi, wann_plot_num), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating wann_func in plot_wannier', comm)
        return
      end if
      wann_func = cmplx_0
      if (spinors) then
        allocate (wann_func_nc(nxx_lo:nxx_hi, nyy_lo:nyy_hi, nzz_lo:nzz_hi, 2, wann_plot_num), &
                  stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating wann_func_nc in plot_wannier', comm)
          return
        end if
        wann_func_nc = cmplx_0
      end if
      if (.not. spinors) then
        if (have_disentangled) then
          allocate (r_wvfn_tmp(ngx*ngy*ngz, maxval(dis_manifold%ndimwin)), stat=ierr)
          if (ierr /= 0) then
            call set_error_alloc(error, 'Error in allocating r_wvfn_tmp in plot_wannier', comm)
            return
          end if
        end if
        allocate (r_wvfn(ngx*ngy*ngz, num_wann), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating r_wvfn in plot_wannier', comm)
          return
        end if
      else
        if (have_disentangled) then
          allocate (r_wvfn_tmp_nc(ngx*ngy*ngz, maxval(dis_manifold%ndimwin), 2), stat=ierr)
          if (ierr /= 0) then
            call set_error_alloc(error, 'Error in allocating r_wvfn_tmp_nc in plot_wannier', comm)
            return
          end if
        end if
        allocate (r_wvfn_nc(ngx*ngy*ngz, num_wann, 2), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating r_wvfn_nc in plot_wannier', comm)
          return
        end if
      end if

      if (.not. spinors) then
        allocate (c_wvfn(ngrid, wann_plot_num), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating c_wvfn in plot_wannier', comm)
          return
        end if
      else
        allocate (c_wvfn_nc(ngrid, wann_plot_num, 2), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating c_wvfn_nc in plot_wannier', comm)
          return
        end if
      end if

      allocate (phase_x(nxx_lo:nxx_hi), phase_y(nyy_lo:nyy_hi), phase_z(nzz_lo:nzz_hi), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating phase arrays in plot_wannier', comm)
        return
      end if

      call io_date(cdate, ctime)
      do loop_kpt = 1, num_kpts
        if (dist_k(loop_kpt) /= my_node_id) cycle

        inc_band = .true.
        num_inc = num_wann
        if (have_disentangled) then
          inc_band(:) = dis_manifold%lwindow(:, loop_kpt)
          num_inc = dis_manifold%ndimwin(loop_kpt)
        end if

        if (.not. spinors) then
          write (wfnname, 200) loop_kpt, wvfn_read%spin_channel
        else
          write (wfnname, 199) loop_kpt
        end if

        if (have_disentangled) then
          call plot_read_unk(wfnname, wvfn_read%formatted, spinors, inc_band, num_inc, &
                             num_bands, loop_kpt, ngx, ngy, ngz, r_wvfn_tmp, r_wvfn_tmp_nc, &
                             stdout, error, comm)
        else
          call plot_read_unk(wfnname, wvfn_read%formatted, spinors, inc_band, num_bands, &
                             num_bands, loop_kpt, ngx, ngy, ngz, r_wvfn, r_wvfn_nc, &
                             stdout, error, comm)
        end if
        if (allocated(error)) return

        if (have_disentangled) then
          if (.not. spinors) then
            r_wvfn = cmplx_0
            do loop_w = 1, num_wann
              do loop_b = 1, num_inc
                r_wvfn(:, loop_w) = r_wvfn(:, loop_w) + &
                                    u_matrix_opt(loop_b, loop_w, loop_kpt)*r_wvfn_tmp(:, loop_b)
              end do
            end do
          else
            r_wvfn_nc = cmplx_0
            do loop_w = 1, num_wann
              do loop_b = 1, num_inc
                call zaxpy(ngx*ngy*ngz, u_matrix_opt(loop_b, loop_w, loop_kpt), r_wvfn_tmp_nc(1, loop_b, 1), 1, & ! up-spinor
                           r_wvfn_nc(1, loop_w, 1), 1)
                call zaxpy(ngx*ngy*ngz, u_matrix_opt(loop_b, loop_w, loop_kpt), r_wvfn_tmp_nc(1, loop_b, 2), 1, & ! down-spinor
                           r_wvfn_nc(1, loop_w, 2), 1)
              end do
            end do
          end if
        end if

        ! Contract band index: c_wvfn(npoint, w) = sum_b u_matrix(b, list(w), k) * r_wvfn(npoint, b)
        if (.not. spinors) then
          c_wvfn = cmplx_0
        else
          c_wvfn_nc = cmplx_0
        end if
        do loop_b = 1, num_wann
          do loop_w = 1, wann_plot_num
            uw = u_matrix(loop_b, wannier_plot%list(loop_w), loop_kpt)
            do npoint = 1, ngrid
              if (.not. spinors) then
                c_wvfn(npoint, loop_w) = c_wvfn(npoint, loop_w) + uw*r_wvfn(npoint, loop_b)
              else
                c_wvfn_nc(npoint, loop_w, 1) = c_wvfn_nc(npoint, loop_w, 1) + uw*r_wvfn_nc(npoint, loop_b, 1)
                c_wvfn_nc(npoint, loop_w, 2) = c_wvfn_nc(npoint, loop_w, 2) + uw*r_wvfn_nc(npoint, loop_b, 2)
              end if
            end do
          end do
        end do

        ! Precompute factored phase arrays for this k-point
        do nxx = nxx_lo, nxx_hi
          phase_x(nxx) = exp(twopi*cmplx_i*kpt_latt(1, loop_kpt)*real(nxx - 1, dp)/real(ngx, dp))
        end do
        do nyy = nyy_lo, nyy_hi
          phase_y(nyy) = exp(twopi*cmplx_i*kpt_latt(2, loop_kpt)*real(nyy - 1, dp)/real(ngy, dp))
        end do
        do nzz = nzz_lo, nzz_hi
          phase_z(nzz) = exp(twopi*cmplx_i*kpt_latt(3, loop_kpt)*real(nzz - 1, dp)/real(ngz, dp))
        end do

        ! nxx, nyy, nzz span a parallelogram in the real space mesh, of side
        ! 2*nphir, and centered around the maximum of phi_i, nphimx(i, 1 2 3)
        !
        ! nx ny nz are the nxx nyy nzz brought back to the unit cell in
        ! which u_nk(r)=cptwrb(r,n)  is represented
        !
        ! There is a big performance improvement in looping over num_wann
        ! in the inner loop. This is poor memory access for wann_func and
        ! but the reduced number of operations wins out.
        do nzz = nzz_lo, nzz_hi
          nz = mod(nzz, ngz)
          if (nz .lt. 1) nz = nz + ngz
          do nyy = nyy_lo, nyy_hi
            ny = mod(nyy, ngy)
            if (ny .lt. 1) ny = ny + ngy
            phase_yz = phase_y(nyy)*phase_z(nzz)
            do nxx = nxx_lo, nxx_hi
              nx = mod(nxx, ngx)
              if (nx .lt. 1) nx = nx + ngx
              npoint = nx + (ny - 1)*ngx + (nz - 1)*ngy*ngx
              catmp = phase_x(nxx)*phase_yz
              do loop_w = 1, wann_plot_num
                if (.not. spinors) then
                  wann_func(nxx, nyy, nzz, loop_w) = &
                    wann_func(nxx, nyy, nzz, loop_w) + c_wvfn(npoint, loop_w)*catmp
                else
                  wann_func_nc(nxx, nyy, nzz, 1, loop_w) = &
                    wann_func_nc(nxx, nyy, nzz, 1, loop_w) + c_wvfn_nc(npoint, loop_w, 1)*catmp
                  wann_func_nc(nxx, nyy, nzz, 2, loop_w) = &
                    wann_func_nc(nxx, nyy, nzz, 2, loop_w) + c_wvfn_nc(npoint, loop_w, 2)*catmp
                end if
              end do
            end do
          end do
        end do

      end do !loop over kpoints

      if (spinors) then
        call comms_reduce(wann_func_nc(nxx_lo, nyy_lo, nzz_lo, 1, 1), &
                          size(wann_func_nc), 'SUM', error, comm)
      else
        call comms_reduce(wann_func(nxx_lo, nyy_lo, nzz_lo, 1), &
                          size(wann_func), 'SUM', error, comm)
      end if
      if (allocated(error)) return

      if (on_root) then
        if (spinors) then
          do nzz = nzz_lo, nzz_hi
            do nyy = nyy_lo, nyy_hi
              do nxx = nxx_lo, nxx_hi
                do loop_w = 1, wann_plot_num
                  upspinor = real(wann_func_nc(nxx, nyy, nzz, 1, loop_w)* &
                                  conjg(wann_func_nc(nxx, nyy, nzz, 1, loop_w)), dp)
                  dnspinor = real(wann_func_nc(nxx, nyy, nzz, 2, loop_w)* &
                                  conjg(wann_func_nc(nxx, nyy, nzz, 2, loop_w)), dp)
                  if (wannier_plot%spinor_phase) then
                    upphase = sign(1.0_dp, real(wann_func_nc(nxx, nyy, nzz, 1, loop_w), dp))
                    dnphase = sign(1.0_dp, real(wann_func_nc(nxx, nyy, nzz, 2, loop_w), dp))
                  else
                    upphase = 1.0_dp; dnphase = 1.0_dp
                  end if
                  select case (wannier_plot%spinor_mode)
                  case ('total')
                    wann_func(nxx, nyy, nzz, loop_w) = cmplx(sqrt(upspinor + dnspinor), 0.0_dp, dp)
                  case ('up')
                    wann_func(nxx, nyy, nzz, loop_w) = cmplx(sqrt(upspinor), 0.0_dp, dp)*upphase
                  case ('down')
                    wann_func(nxx, nyy, nzz, loop_w) = cmplx(sqrt(dnspinor), 0.0_dp, dp)*dnphase
                  case default
                    call set_error_file(error, 'plot_wannier: Invalid wannier_plot_spinor_mode ' &
                                        //trim(wannier_plot%spinor_mode), comm)
                    return
                  end select
                  wann_func(nxx, nyy, nzz, loop_w) = wann_func(nxx, nyy, nzz, loop_w)/real(num_kpts, dp)
                end do
              end do
            end do
          end do
        end if
      end if

      if (on_root) then
        if (.not. spinors) then !!!!! For spinor Wannier functions, the steps below are not necessary.
          ! fix the global phase by setting the wannier to
          ! be real at the point where it has max. modulus

          do loop_w = 1, wann_plot_num
            tmaxx = 0.0
            wmod = cmplx_1
            do nzz = nzz_lo, nzz_hi
              do nyy = nyy_lo, nyy_hi
                do nxx = nxx_lo, nxx_hi
                  wann_func(nxx, nyy, nzz, loop_w) = wann_func(nxx, nyy, nzz, loop_w)/real(num_kpts, dp)
                  tmax = real(wann_func(nxx, nyy, nzz, loop_w)* &
                              conjg(wann_func(nxx, nyy, nzz, loop_w)), dp)
                  if (tmax > tmaxx) then
                    tmaxx = tmax
                    wmod = wann_func(nxx, nyy, nzz, loop_w)
                  end if
                end do
              end do
            end do
            wmod = wmod/sqrt(real(wmod)**2 + aimag(wmod)**2)
            wann_func(:, :, :, loop_w) = wann_func(:, :, :, loop_w)/wmod
            write (stdout, '(6x,a,i4,7x,a,f11.6,SP,f11.6,"i")') "Wannier Function Num: ", &
              wannier_plot%list(loop_w), "Phase Factor = ", 1/wmod
          end do
          write (stdout, *) ''
          !
          ! Check the 'reality' of the WF
          !
          do loop_w = 1, wann_plot_num
            ratmax = 0.0_dp
            do nzz = nzz_lo, nzz_hi
              do nyy = nyy_lo, nyy_hi
                do nxx = nxx_lo, nxx_hi
                  if (abs(real(wann_func(nxx, nyy, nzz, loop_w), dp)) >= 0.01_dp) then
                    ratio = abs(aimag(wann_func(nxx, nyy, nzz, loop_w)))/ &
                            abs(real(wann_func(nxx, nyy, nzz, loop_w), dp))
                    ratmax = max(ratmax, ratio)
                  end if
                end do
              end do
            end do
            write (stdout, '(6x,a,i4,7x,a,f11.6)') 'Wannier Function Num: ', wannier_plot%list(loop_w), &
              'Maximum Im/Re Ratio = ', ratmax
          end do
        end if !!!!!
        write (stdout, *) ' '
        if (wannier_plot%format .eq. 'xcrysden') then
          call internal_xsf_format()
        elseif (wannier_plot%format .eq. 'cube') then
          call internal_cube_format(atom_data, wannier_data, real_lattice, bohr, error)
          if (allocated(error)) return
        else
          call set_error_warn(error, 'wannier_plot_format not recognised in plot_wannier', comm)
          return
        end if

        if (print_output%timing_level > 1) call io_stopwatch_stop('plot: wannier', timer)
      end if !on_root

    end associate

    if (allocated(c_wvfn)) then
      deallocate (c_wvfn, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating c_wvfn in plot_wannier', comm)
        return
      end if
    end if
    if (allocated(c_wvfn_nc)) then
      deallocate (c_wvfn_nc, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating c_wvfn_nc in plot_wannier', comm)
        return
      end if
    end if
    if (allocated(phase_x)) then
      deallocate (phase_x, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating phase_x in plot_wannier', comm)
        return
      end if
    end if
    if (allocated(phase_y)) then
      deallocate (phase_y, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating phase_y in plot_wannier', comm)
        return
      end if
    end if
    if (allocated(phase_z)) then
      deallocate (phase_z, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating phase_z in plot_wannier', comm)
        return
      end if
    end if
    return

  contains

    !================================================!
    subroutine internal_cube_format(atom_data, wannier_data, real_lattice, bohr, error)
      !================================================!
      !
      !! Write WFs in Gaussian cube format.
      !
      !================================================!

      use w90_utility, only: utility_translate_home, utility_cart_to_frac, utility_frac_to_cart, &
                             utility_inverse_mat, utility_recip_lattice_base
      use w90_types, only: wannier_data_type, atom_data_type
      use w90_wannier90_types, only: wvfn_read_type
      use w90_error, only: w90_error_type, set_error_alloc, set_error_warn, set_error_dealloc

      implicit none

      type(wannier_data_type), intent(in) :: wannier_data
      type(atom_data_type), intent(in) :: atom_data
      type(w90_error_type), allocatable, intent(out) :: error
      real(kind=dp), intent(in) :: bohr

      real(kind=dp), intent(in) :: real_lattice(3, 3)

      real(kind=dp), allocatable :: wann_cube(:, :, :)
      real(kind=dp) :: inv_lattice(3, 3), recip_lattice(3, 3), pos_frac(3), volume
      real(kind=dp) :: rstart(3), rend(3), rlength(3), orig(3), dgrid(3)
      real(kind=dp) :: moda(3), modb(3)
      real(kind=dp) :: val_Q
      real(kind=dp) :: comf(3), wcf(3), diff(3), difc(3), dist
      integer :: ierr, iname, max_elements
      integer :: isp, iat, nzz, nyy, nxx, loop_w, qxx, qyy, qzz, wann_index
      integer :: istart(3), iend(3), ilength(3), nend
      integer :: ixx, iyy, izz
      integer :: irdiff(3), icount
      integer, allocatable :: atomic_Z(:)
      logical :: lmol, lcrys
      character(len=2), dimension(109) :: periodic_table = (/ &
           & 'h ', 'he', &
           & 'li', 'be', 'b ', 'c ', 'n ', 'o ', 'f ', 'ne', &
           & 'na', 'mg', 'al', 'si', 'p ', 's ', 'cl', 'ar', &
           & 'k ', 'ca', 'sc', 'ti', 'v ', 'cr', 'mn', 'fe', 'co', 'ni', 'cu', 'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr', &
           & 'rb', 'sr', 'y ', 'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag', 'cd', 'in', 'sn', 'sb', 'te', 'i ', 'xe', &
           & 'cs', 'ba', &
           & 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy', 'ho', 'er', 'tm', 'yb', 'lu', &
           & 'hf', 'ta', 'w ', 're', 'os', 'ir', 'pt', 'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn', &
           & 'fr', 'ra', &
           & 'ac', 'th', 'pa', 'u ', 'np', 'pu', 'am', 'cm', 'bk', 'cf', 'es', 'fm', 'md', 'no', 'lr', &
           & 'rf', 'db', 'sg', 'bh', 'hs', 'mt'/)

      associate (ngs => wannier_plot%supercell)

        allocate (atomic_Z(atom_data%num_species), stat=ierr)
        if (ierr .ne. 0) then
          call set_error_alloc(error, 'Error: allocating atomic_Z in plot_wannier: cube_format', comm)
          return
        end if

        call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
        lmol = .false.
        lcrys = .false.
        if (index(wannier_plot%mode, 'mol') > 0) lmol = .true.      ! molecule mode
        if (index(wannier_plot%mode, 'crys') > 0) lcrys = .true.    ! crystal mode

        val_Q = 1.0_dp ! dummy value for cube file

        ! Assign atomic numbers to species
        max_elements = size(periodic_table)
        atomic_Z(:) = 0
        do isp = 1, atom_data%num_species
          do iname = 1, max_elements
            if (atom_data%symbol(isp) .eq. periodic_table(iname)) then
              atomic_Z(isp) = iname
              exit
            end if
          end do
        end do

202     format(a, '_', i5.5, '.cube')

        ! Lengths of real and reciprocal lattice vectors
        do i = 1, 3
          moda(i) = sqrt(real_lattice(i, 1)*real_lattice(i, 1) &
                         + real_lattice(i, 2)*real_lattice(i, 2) &
                         + real_lattice(i, 3)*real_lattice(i, 3))
          modb(i) = sqrt(recip_lattice(i, 1)*recip_lattice(i, 1) &
                         + recip_lattice(i, 2)*recip_lattice(i, 2) &
                         + recip_lattice(i, 3)*recip_lattice(i, 3))
        end do

        ! Grid spacing in each lattice direction
        dgrid(1) = moda(1)/ngx; dgrid(2) = moda(2)/ngy; dgrid(3) = moda(3)/ngz

        ! Find "centre of mass" of atomic positions (in fractional coordinates)
        call utility_inverse_mat(real_lattice, inv_lattice)
        comf(:) = 0.0_dp
        do isp = 1, atom_data%num_species
          do iat = 1, atom_data%species_num(isp)
            call utility_cart_to_frac(atom_data%pos_cart(:, iat, isp), pos_frac, inv_lattice)
            comf(:) = comf(:) + pos_frac(:)
          end do
        end do
        comf(:) = comf(:)/atom_data%num_atoms

        ! Loop over WFs
        do loop_w = 1, wann_plot_num

          wann_index = wannier_plot%list(loop_w)
          write (wancube, 202) trim(seedname), wann_index

          ! Find start and end of cube wrt simulation (home) cell origin
          do i = 1, 3
            ! ... in terms of distance along each lattice vector direction i
            rstart(i) = (wannier_data%centres(1, wann_index)*recip_lattice(i, 1) &
                         + wannier_data%centres(2, wann_index)*recip_lattice(i, 2) &
                         + wannier_data%centres(3, wann_index)*recip_lattice(i, 3))*moda(i)/twopi &
                        - twopi*wannier_plot%radius/(moda(i)*modb(i))
            rend(i) = (wannier_data%centres(1, wann_index)*recip_lattice(i, 1) &
                       + wannier_data%centres(2, wann_index)*recip_lattice(i, 2) &
                       + wannier_data%centres(3, wann_index)*recip_lattice(i, 3))*moda(i)/twopi &
                      + twopi*wannier_plot%radius/(moda(i)*modb(i))
          end do

          rlength(:) = rend(:) - rstart(:)
          ilength(:) = ceiling(rlength(:)/dgrid(:))

          ! ... in terms of integer gridpoints along each lattice vector direction i
          istart(:) = floor(rstart(:)/dgrid(:)) + 1
          iend(:) = istart(:) + ilength(:) - 1

          ! Origin of cube wrt simulation (home) cell in Cartesian co-ordinates
          do i = 1, 3
            orig(i) = real(istart(1) - 1, dp)*dgrid(1)*real_lattice(1, i)/moda(1) &
                      + real(istart(2) - 1, dp)*dgrid(2)*real_lattice(2, i)/moda(2) &
                      + real(istart(3) - 1, dp)*dgrid(3)*real_lattice(3, i)/moda(3)
          end do

          ! Debugging
          if (print_output%iprint > 3) then
            write (stdout, '(a,i12)') 'loop_w  =', loop_w
            write (stdout, '(a,3f12.6)') 'comf    =', (comf(i), i=1, 3)
            write (stdout, '(a,3i12)') 'ngi     =', ngx, ngy, ngz
            write (stdout, '(a,3f12.6)') 'dgrid   =', (dgrid(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'rstart  =', (rstart(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'rend    =', (rend(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'rlength =', (rlength(i), i=1, 3)
            write (stdout, '(a,3i12)') 'istart  =', (istart(i), i=1, 3)
            write (stdout, '(a,3i12)') 'iend    =', (iend(i), i=1, 3)
            write (stdout, '(a,3i12)') 'ilength =', (ilength(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'orig    =', (orig(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'wann_cen=', (wannier_data%centres(i, wann_index), i=1, 3)
          end if

          allocate (wann_cube(1:ilength(1), 1:ilength(2), 1:ilength(3)), stat=ierr)
          if (ierr .ne. 0) then
            call set_error_alloc(error, 'Error: allocating wann_cube in plot_wannier: cube_format', comm)
            return
          end if

          ! initialise
          wann_cube = 0.0_dp

          do nzz = 1, ilength(3)
            qzz = nzz + istart(3) - 1
            izz = int((abs(qzz) - 1)/ngz)
!            if (qzz.lt.-ngz) qzz=qzz+izz*ngz
!            if (qzz.gt.(ngs(3)-1)*ngz-1) then
            if (qzz .lt. nzz_lo) qzz = qzz + izz*ngz
            if (qzz .gt. nzz_hi) then
              write (stdout, *) 'Error plotting WF cube. Try one of the following:'
              write (stdout, *) '   (1) increase wannier_plot_supercell;'
              write (stdout, *) '   (2) decrease wannier_plot_radius;'
              write (stdout, *) '   (3) set wannier_plot_format=xcrysden'
              call set_error_warn(error, 'Error plotting WF cube.', comm)
              return
            end if
            do nyy = 1, ilength(2)
              qyy = nyy + istart(2) - 1
              iyy = int((abs(qyy) - 1)/ngy)
!               if (qyy.lt.-ngy) qyy=qyy+iyy*ngy
!               if (qyy.gt.(ngs(2)-1)*ngy-1) then
              if (qyy .lt. nyy_lo) qyy = qyy + iyy*ngy
              if (qyy .gt. nyy_hi) then
                write (stdout, *) 'Error plotting WF cube. Try one of the following:'
                write (stdout, *) '   (1) increase wannier_plot_supercell;'
                write (stdout, *) '   (2) decrease wannier_plot_radius;'
                write (stdout, *) '   (3) set wannier_plot_format=xcrysden'
                call set_error_warn(error, 'Error plotting WF cube.', comm)
                return
              end if
              do nxx = 1, ilength(1)
                qxx = nxx + istart(1) - 1
                ixx = int((abs(qxx) - 1)/ngx)
!                  if (qxx.lt.-ngx) qxx=qxx+ixx*ngx
!                  if (qxx.gt.(ngs(1)-1)*ngx-1) then
                if (qxx .lt. nxx_lo) qxx = qxx + ixx*ngx
                if (qxx .gt. nxx_hi) then
                  write (stdout, *) 'Error plotting WF cube. Try one of the following:'
                  write (stdout, *) '   (1) increase wannier_plot_supercell;'
                  write (stdout, *) '   (2) decrease wannier_plot_radius;'
                  write (stdout, *) '   (3) set wannier_plot_format=xcrysden'
                  call set_error_warn(error, 'Error plotting WF cube.', comm)
                  return
                end if
                wann_cube(nxx, nyy, nzz) = real(wann_func(qxx, qyy, qzz, loop_w), dp)
              end do
            end do
          end do

          ! WF centre in fractional coordinates
          call utility_cart_to_frac(wannier_data%centres(:, wann_index), wcf(:), inv_lattice)

          ! The vector (in fractional coordinates) from WF centre to "centre of mass"
          diff(:) = comf(:) - wcf(:)

          ! Corresponding nearest cell vector
          irdiff(:) = nint(diff(:))

          if (print_output%iprint > 3) then
            write (stdout, '(a,3f12.6)') 'wcf     =', (wcf(i), i=1, 3)
            write (stdout, '(a,3f12.6)') 'diff    =', (diff(i), i=1, 3)
            write (stdout, '(a,3i12)') 'irdiff  =', (irdiff(i), i=1, 3)
          end if

          if (lmol) then ! In "molecule mode" translate origin of cube to bring it in coincidence with the atomic positions
            orig(:) = orig(:) + real(irdiff(1), kind=dp)*real_lattice(1, :) &
                      + real(irdiff(2), kind=dp)*real_lattice(2, :) &
                      + real(irdiff(3), kind=dp)*real_lattice(3, :)
            if (print_output%iprint > 3) write (stdout, '(a,3f12.6,/)') 'orig-new=', (orig(i), i=1, 3)
          else ! In "crystal mode" count number of atoms within a given radius of wannier centre
            icount = 0
            do isp = 1, atom_data%num_species
              do iat = 1, atom_data%species_num(isp)
                call utility_cart_to_frac(atom_data%pos_cart(:, iat, isp), pos_frac, inv_lattice)
                do nzz = -ngs(3)/2, (ngs(3) + 1)/2
                  do nyy = -ngs(2)/2, (ngs(2) + 1)/2
                    do nxx = -ngs(1)/2, (ngs(1) + 1)/2
                      diff(:) = pos_frac(:) - wcf(:) &
                                + (/real(nxx, kind=dp), real(nyy, kind=dp), real(nzz, kind=dp)/)
                      call utility_frac_to_cart(diff, difc, real_lattice)
                      dist = sqrt(difc(1)*difc(1) + difc(2)*difc(2) + difc(3)*difc(3))
                      if (dist .le. (wannier_plot%scale*wannier_plot%radius)) then
                        icount = icount + 1
                      end if
                    end do
                  end do
                end do
              end do ! iat
            end do ! isp
            if (print_output%iprint > 3) write (stdout, '(a,i12)') 'icount  =', icount
          end if

          ! Write cube file (everything in Bohr)
          open (newunit=file_unit, file=trim(wancube), form='formatted', status='unknown')
          ! First two lines are comments
          write (file_unit, *) '     Generated by Wannier90 code http://www.wannier.org'
          write (file_unit, *) '     On ', cdate, ' at ', ctime
          ! Number of atoms, origin of cube (Cartesians) wrt simulation (home) cell
          if (lmol) then
            write (file_unit, '(i4,3f13.5)') atom_data%num_atoms, orig(1)/bohr, orig(2)/bohr, orig(3)/bohr
          else
            write (file_unit, '(i4,3f13.5)') icount, orig(1)/bohr, orig(2)/bohr, orig(3)/bohr
          end if
          ! Number of grid points in each direction, lattice vector
          write (file_unit, '(i4,3f13.5)') ilength(1), real_lattice(1, 1)/(real(ngx, dp)*bohr), &
            real_lattice(1, 2)/(real(ngx, dp)*bohr), real_lattice(1, 3)/(real(ngx, dp)*bohr)
          write (file_unit, '(i4,3f13.5)') ilength(2), real_lattice(2, 1)/(real(ngy, dp)*bohr), &
            real_lattice(2, 2)/(real(ngy, dp)*bohr), real_lattice(2, 3)/(real(ngy, dp)*bohr)
          write (file_unit, '(i4,3f13.5)') ilength(3), real_lattice(3, 1)/(real(ngz, dp)*bohr), &
            real_lattice(3, 2)/(real(ngz, dp)*bohr), real_lattice(3, 3)/(real(ngz, dp)*bohr)

          ! Atomic number, valence charge, position of atom
!         do isp=1,num_species
!            do iat=1,atoms_species_num(isp)
!               write(file_unit,'(i4,4f13.5)') atomic_Z(isp), val_Q, (atoms_pos_cart(i,iat,isp)/bohr,i=1,3)
!            end do
!         end do

          do isp = 1, atom_data%num_species
            do iat = 1, atom_data%species_num(isp)
              if (lmol) then ! In "molecule mode", write atomic coordinates as they appear in input file
                write (file_unit, '(i4,4f13.5)') atomic_Z(isp), val_Q, (atom_data%pos_cart(i, iat, isp)/bohr, i=1, 3)
              else           ! In "crystal mode", write atoms in supercell within a given radius of Wannier centre
                call utility_cart_to_frac(atom_data%pos_cart(:, iat, isp), pos_frac, inv_lattice)
                do nzz = -ngs(3)/2, (ngs(3) + 1)/2
                  do nyy = -ngs(2)/2, (ngs(2) + 1)/2
                    do nxx = -ngs(1)/2, (ngs(1) + 1)/2
                      diff(:) = pos_frac(:) - wcf(:) &
                                + (/real(nxx, kind=dp), real(nyy, kind=dp), real(nzz, kind=dp)/)
                      call utility_frac_to_cart(diff, difc, real_lattice)
                      dist = sqrt(difc(1)*difc(1) + difc(2)*difc(2) + difc(3)*difc(3))
                      if (dist .le. (wannier_plot%scale*wannier_plot%radius)) then
                        diff(:) = pos_frac(:) &
                                  + (/real(nxx, kind=dp), real(nyy, kind=dp), real(nzz, kind=dp)/)
                        call utility_frac_to_cart(diff, difc, real_lattice)
                        write (file_unit, '(i4,4f13.5)') atomic_Z(isp), val_Q, (difc(i)/bohr, i=1, 3)
                      end if
                    end do
                  end do
                end do
              end if
            end do ! iat
          end do ! isp

          ! Volumetric data in batches of 6 values per line, 'z'-direction first.
          do nxx = 1, ilength(1)
            do nyy = 1, ilength(2)
              do nzz = 1, ilength(3), 6
                nend = min(nzz + 5, ilength(3))
                write (file_unit, '(6E13.5)') wann_cube(nxx, nyy, nzz:nend)
              end do
            end do
          end do

          deallocate (wann_cube, stat=ierr)
          if (ierr .ne. 0) then
            call set_error_dealloc(error, 'Error: deallocating wann_cube in plot_wannier: cube_format', comm)
            return
          end if

        end do

        deallocate (atomic_Z, stat=ierr)
        if (ierr .ne. 0) then
          call set_error_dealloc(error, 'Error: deallocating atomic_Z in plot_wannier: cube_format', comm)
          return
        end if

      end associate

      return

    end subroutine internal_cube_format

    subroutine internal_xsf_format()

      implicit none

201   format(a, '_', i5.5, '.xsf')

      associate (ngs => wannier_plot%supercell)

        ! this is to create the WF...xsf output, to be read by XCrySDen
        ! (coordinates + isosurfaces)

        x_0ang = -real(((ngs(1))/2)*ngx + 1, dp)/real(ngx, dp)*real_lattice(1, 1) - &
                 real(((ngs(2))/2)*ngy + 1, dp)/real(ngy, dp)*real_lattice(2, 1) - &
                 real(((ngs(3))/2)*ngz + 1, dp)/real(ngz, dp)*real_lattice(3, 1)
        y_0ang = -real(((ngs(1))/2)*ngx + 1, dp)/real(ngx, dp)*real_lattice(1, 2) - &
                 real(((ngs(2))/2)*ngy + 1, dp)/real(ngy, dp)*real_lattice(2, 2) - &
                 real(((ngs(3))/2)*ngz + 1, dp)/real(ngz, dp)*real_lattice(3, 2)
        z_0ang = -real(((ngs(1))/2)*ngx + 1, dp)/real(ngx, dp)*real_lattice(1, 3) - &
                 real(((ngs(2))/2)*ngy + 1, dp)/real(ngy, dp)*real_lattice(2, 3) - &
                 real(((ngs(3))/2)*ngz + 1, dp)/real(ngz, dp)*real_lattice(3, 3)

        fxcry(1) = real(ngs(1)*ngx - 1, dp)/real(ngx, dp)
        fxcry(2) = real(ngs(2)*ngy - 1, dp)/real(ngy, dp)
        fxcry(3) = real(ngs(3)*ngz - 1, dp)/real(ngz, dp)
        do j = 1, 3
          dirl(:, j) = fxcry(:)*real_lattice(:, j)
        end do

        do loop_b = 1, wann_plot_num

          write (wanxsf, 201) trim(seedname), wannier_plot%list(loop_b)

          open (newunit=file_unit, file=trim(wanxsf), form='formatted', status='unknown')
          write (file_unit, *) '      #'
          write (file_unit, *) '      # Generated by the Wannier90 code http://www.wannier.org'
          write (file_unit, *) '      # On ', cdate, ' at ', ctime
          write (file_unit, *) '      #'
          ! should pass this into the code
          if (index(wannier_plot%mode, 'mol') > 0) then
            write (file_unit, '("ATOMS")')
          else
            write (file_unit, '("CRYSTAL")')
            write (file_unit, '("PRIMVEC")')
            write (file_unit, '(3f12.7)') real_lattice(1, 1), real_lattice(1, 2), real_lattice(1, 3)
            write (file_unit, '(3f12.7)') real_lattice(2, 1), real_lattice(2, 2), real_lattice(2, 3)
            write (file_unit, '(3f12.7)') real_lattice(3, 1), real_lattice(3, 2), real_lattice(3, 3)
            write (file_unit, '("CONVVEC")')
            write (file_unit, '(3f12.7)') real_lattice(1, 1), real_lattice(1, 2), real_lattice(1, 3)
            write (file_unit, '(3f12.7)') real_lattice(2, 1), real_lattice(2, 2), real_lattice(2, 3)
            write (file_unit, '(3f12.7)') real_lattice(3, 1), real_lattice(3, 2), real_lattice(3, 3)
            write (file_unit, '("PRIMCOORD")')
            write (file_unit, '(i6,"  1")') atom_data%num_atoms
          end if
          do nsp = 1, atom_data%num_species
            do nat = 1, atom_data%species_num(nsp)
              write (file_unit, '(a2,3x,3f12.7)') atom_data%label(nsp), (atom_data%pos_cart(i, nat, nsp), i=1, 3)
            end do
          end do

          write (file_unit, '(/)')
          write (file_unit, '("BEGIN_BLOCK_DATAGRID_3D",/,"3D_field",/, "BEGIN_DATAGRID_3D_UNKNOWN")')
          write (file_unit, '(3i6)') ngs(1)*ngx, ngs(2)*ngy, ngs(3)*ngz
          write (file_unit, '(3f12.6)') x_0ang, y_0ang, z_0ang
          write (file_unit, '(3f12.7)') dirl(1, 1), dirl(1, 2), dirl(1, 3)
          write (file_unit, '(3f12.7)') dirl(2, 1), dirl(2, 2), dirl(2, 3)
          write (file_unit, '(3f12.7)') dirl(3, 1), dirl(3, 2), dirl(3, 3)
          write (file_unit, '(6e13.5)') &
            (((real(wann_func(nx, ny, nz, loop_b)), nx=nxx_lo, nxx_hi), &
              ny=nyy_lo, nyy_hi), nz=nzz_lo, nzz_hi)
          write (file_unit, '("END_DATAGRID_3D",/, "END_BLOCK_DATAGRID_3D")')
          close (file_unit)

        end do

      end associate
      return
    end subroutine internal_xsf_format

  end subroutine plot_wannier