plot_interpolate_bands

private subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, real_space_ham, ws_region, print_output, recip_lattice, num_wann, wannier_data, ham_r, irvec, ndegen, nrpts, wannier_centres_translated, ws_distance, bands_num_spec_points, stdout, seedname, timer, error, comm)

Uses

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Plots the interpolated band structure
Arguments

Type	Intent	Attributes		Name
integer,	intent(in)		::	mp_grid(3)
real(kind=dp),	intent(in)		::	real_lattice(3,3)
type(band_plot_type),	intent(in)		::	band_plot
type(kpoint_path_type),	intent(in)		::	kpoint_path
type(real_space_ham_type),	intent(in)		::	real_space_ham
type(ws_region_type),	intent(in)		::	ws_region
type(print_output_type),	intent(in)		::	print_output
real(kind=dp),	intent(in)		::	recip_lattice(3,3)
integer,	intent(in)		::	num_wann
type(wannier_data_type),	intent(in)		::	wannier_data
complex(kind=dp),	intent(in)		::	ham_r(:,:,:)
integer,	intent(in)		::	irvec(:,:)
integer,	intent(in)		::	ndegen(:)
integer,	intent(inout)		::	nrpts
real(kind=dp),	intent(in)		::	wannier_centres_translated(:,:)
type(ws_distance_type),	intent(inout)		::	ws_distance
integer,	intent(in)		::	bands_num_spec_points
integer,	intent(in)		::	stdout
character(len=50),	intent(in)		::	seedname
type(timer_list_type),	intent(inout)		::	timer
type(w90_error_type),	intent(out),	allocatable	::	error
type(w90_comm_type),	intent(in)		::	comm
Called by

Help
Source Code

  subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, real_space_ham, &
                                    ws_region, print_output, recip_lattice, num_wann, &
                                    wannier_data, ham_r, irvec, ndegen, nrpts, &
                                    wannier_centres_translated, ws_distance, &
                                    bands_num_spec_points, stdout, seedname, timer, error, comm)
    !================================================!
    !                                            !
    !! Plots the interpolated band structure
    !                                            !
    !================================================!

    use w90_constants, only: dp, cmplx_0, twopi
    use w90_io, only: io_time, io_stopwatch_start, io_stopwatch_stop
    use w90_ws_distance, only: ws_translate_dist
    use w90_utility, only: utility_metric
    use w90_types, only: wannier_data_type, kpoint_path_type, print_output_type, ws_region_type, &
                         ws_distance_type, timer_list_type
    use w90_wannier90_types, only: band_plot_type, real_space_ham_type
    use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal, &
                         set_error_warn

    implicit none

    ! arguments
    type(band_plot_type), intent(in) :: band_plot
    type(kpoint_path_type), intent(in) :: kpoint_path
    type(print_output_type), intent(in) :: print_output
    type(real_space_ham_type), intent(in) :: real_space_ham
    type(wannier_data_type), intent(in) :: wannier_data
    type(ws_distance_type), intent(inout) :: ws_distance
    type(ws_region_type), intent(in) :: ws_region
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(inout) :: nrpts
    integer, intent(in) :: ndegen(:)
    integer, intent(in) :: irvec(:, :)
    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_wann
    integer, intent(in) :: bands_num_spec_points
    integer, intent(in) :: stdout

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

    complex(kind=dp), intent(in) :: ham_r(:, :, :)
    character(len=50), intent(in)  :: seedname

    ! local variables
    integer, allocatable :: irvec_cut(:, :)
    integer              :: irvec_max(3)
    integer              :: nrpts_cut
    integer, allocatable :: iwork(:), ifail(:)
    integer              :: info, i, j
    integer              :: irpt, nfound, loop_kpt, counter
    integer              :: loop_spts, total_pts, loop_i, nkp, ideg
    integer              :: num_paths, num_spts, ierr
    integer              :: bndunit, gnuunit, loop_w, loop_p
    integer, allocatable :: kpath_pts(:)
    integer, allocatable :: idx_special_points(:)
    integer, allocatable :: label_idx_special_points(:)

    real(kind=dp)              :: rdotk, vec(3), emin, emax, time0
    real(kind=dp), allocatable :: kpath_len(:)
    real(kind=dp), allocatable :: rwork(:)
    real(kind=dp), allocatable :: xval(:)
    real(kind=dp), allocatable :: eig_int(:, :), plot_kpoint(:, :)
    real(kind=dp), allocatable :: bands_proj(:, :)
    real(kind=dp), allocatable :: xval_special_points(:)
    real(kind=dp)              :: recip_metric(3, 3)

    complex(kind=dp)              :: fac
    complex(kind=dp), allocatable :: ham_r_cut(:, :, :)
    complex(kind=dp), allocatable :: ham_pack(:)
    complex(kind=dp), allocatable :: ham_kprm(:, :)
    complex(kind=dp), allocatable :: U_int(:, :)
    complex(kind=dp), allocatable :: cwork(:)

    logical, allocatable :: kpath_print_first_point(:)

    character(len=20), allocatable :: glabel(:)
    character(len=20), allocatable :: xlabel(:)
    character(len=20), allocatable :: ctemp(:)

    ! mpi variables
    integer ::  my_node_id, num_nodes
    logical ::  on_root
    integer, allocatable :: counts(:)
    integer, allocatable :: displs(:)

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

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

    allocate (counts(0:num_nodes - 1), stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating counts in plot_interpolate_bands', comm)
      return
    end if
    allocate (displs(0:num_nodes - 1), stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating displs in plot_interpolate_bands', comm)
      return
    end if
    !
    if (on_root) then
      if (print_output%timing_level > 1) then
        call io_stopwatch_start('plot: interpolate_bands', timer)
      end if
      !
      time0 = io_time()

      write (stdout, *)
      write (stdout, '(1x,a)') 'Calculating interpolated band-structure'
      write (stdout, *)
    end if ! on_root
    call utility_metric(recip_lattice, recip_metric)
    !
    allocate (ham_pack((num_wann*(num_wann + 1))/2), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ham_pack in plot_interpolate_bands', comm)
      return
    end if
    allocate (ham_kprm(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ham_kprm in plot_interpolate_bands', comm)
      return
    end if
    allocate (U_int(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating U_int in plot_interpolate_bands', comm)
      return
    end if
    allocate (cwork(2*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cwork in plot_interpolate_bands', comm)
      return
    end if
    allocate (rwork(7*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating rwork in plot_interpolate_bands', comm)
      return
    end if
    allocate (iwork(5*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating iwork in plot_interpolate_bands', comm)
      return
    end if
    allocate (ifail(num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ifail in plot_interpolate_bands', comm)
      return
    end if

    if (.not. kpoint_path%bands_kpt_explicit) then
      allocate (kpath_len(bands_num_spec_points/2), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kpath_len in plot_interpolate_bands', comm)
        return
      end if
      allocate (kpath_pts(bands_num_spec_points/2), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kpath_pts in plot_interpolate_bands', comm)
        return
      end if
      allocate (kpath_print_first_point(bands_num_spec_points/2), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kpath_print_first_point in plot_interpolate_bands', comm)
        return
      end if
      allocate (idx_special_points(bands_num_spec_points), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating idx_special_points in plot_interpolate_bands', comm)
        return
      end if
      allocate (xval_special_points(bands_num_spec_points), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating xval_special_points in plot_interpolate_bands', comm)
        return
      end if
      idx_special_points = -1
      xval_special_points = -1._dp
    end if
    !
    ! Work out how many points in the total path and the positions of the special points
    !
    if (kpoint_path%bands_kpt_explicit) then
      total_pts = size(kpoint_path%bands_kpt_frac, 2)
      ! Count the total number of special points
      num_spts = 0
      do i = 1, total_pts
        do j = 1, bands_num_spec_points
          if (sum((kpoint_path%bands_kpt_frac(:, i) - kpoint_path%points(:, j))**2) <= 1.e-6) then
            num_spts = num_spts + 1
            exit
          end if
        end do
      end do
    else
      num_paths = bands_num_spec_points/2

      kpath_print_first_point = .false.

      ! Loop over paths, set to False print_first_point if the starting point
      ! is the same as the ending point of the previous path.
      ! I skip the first path for which I always want to print the first point.
      kpath_print_first_point(1) = .true.
      do i = 2, num_paths
        ! If either the coordinates are different or the label is different, compute again the point
        ! (it will end up at the same x coordinate)
        if ((SUM((kpoint_path%points(:, (i - 1)*2) - &
                  kpoint_path%points(:, (i - 1)*2 + 1))**2) > 1.e-6) .or. &
            (TRIM(kpoint_path%labels((i - 1)*2)) .ne. &
             TRIM(kpoint_path%labels((i - 1)*2 + 1)))) then
          kpath_print_first_point(i) = .true.
        end if
      end do

      ! Count the total number of special points
      num_spts = num_paths
      do i = 1, num_paths
        if (kpath_print_first_point(i)) num_spts = num_spts + 1
      end do

      do loop_spts = 1, num_paths
        vec = kpoint_path%points(:, 2*loop_spts) - kpoint_path%points(:, 2*loop_spts - 1)
        kpath_len(loop_spts) = sqrt(dot_product(vec, (matmul(recip_metric, vec))))
        if (loop_spts == 1) then
          kpath_pts(loop_spts) = kpoint_path%num_points_first_segment
        else
          kpath_pts(loop_spts) = nint(real(kpoint_path%num_points_first_segment, dp) &
                                      *kpath_len(loop_spts)/kpath_len(1))
          ! At least 1 point
          !if (kpath_pts(loop_spts) .eq. 0) kpath_pts(loop_spts) = 1
        end if
      end do
      total_pts = sum(kpath_pts)
      do i = 1, num_paths
        if (kpath_print_first_point(i)) total_pts = total_pts + 1
      end do
    end if

    allocate (plot_kpoint(3, total_pts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating plot_kpoint in plot_interpolate_bands', comm)
      return
    end if
    allocate (xval(total_pts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating xval in plot_interpolate_bands', comm)
      return
    end if
    allocate (eig_int(num_wann, total_pts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating eig_int in plot_interpolate_bands', comm)
      return
    end if
    allocate (bands_proj(num_wann, total_pts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating bands_proj in plot_interpolate_bands', comm)
      return
    end if
    allocate (glabel(num_spts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating glabel in plot_interpolate_bands', comm)
      return
    end if
    allocate (xlabel(num_spts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating xlabel in plot_interpolate_bands', comm)
      return
    end if
    allocate (ctemp(bands_num_spec_points), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ctemp in plot_interpolate_bands', comm)
      return
    end if
    eig_int = 0.0_dp; bands_proj = 0.0_dp
    !
    ! Find the position of each kpoint in the path
    !
    if (kpoint_path%bands_kpt_explicit) then
      allocate (idx_special_points(num_spts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating idx_special_points in plot_interpolate_bands', comm)
        return
      end if
      allocate (xval_special_points(num_spts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating xval_special_points in plot_interpolate_bands', comm)
        return
      end if
      idx_special_points = -1
      xval_special_points = -1._dp
      allocate (label_idx_special_points(num_spts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating label_idx_special_points in plot_interpolate_bands', comm)
        return
      end if
      plot_kpoint(:, :) = kpoint_path%bands_kpt_frac(:, :)
      xval = 0.0_dp
      counter = 0
      do i = 1, total_pts
        if (i > 1) then
          vec = plot_kpoint(:, i) - plot_kpoint(:, i - 1)
          xval(i) = xval(i - 1) + sqrt(dot_product(vec, (matmul(recip_metric, vec))))
        end if
        do j = 1, bands_num_spec_points
          if (sum((kpoint_path%bands_kpt_frac(:, i) - kpoint_path%points(:, j))**2) <= 1.e-6) then
            counter = counter + 1
            idx_special_points(counter) = i
            label_idx_special_points(counter) = j
            xval_special_points(counter) = xval(i)
            if (counter > 1) then
              if (idx_special_points(counter) == idx_special_points(counter - 1) + 1) then
                ! If the two points are consecutive, it means that the x coordinate should be the same
                xval(i) = xval(i - 1)
                xval_special_points(counter) = xval(i)
              end if
            end if
            exit
          end if
        end do
      end do

    else
      counter = 0
      do loop_spts = 1, num_paths
        if (kpath_print_first_point(loop_spts)) then
          counter = counter + 1
          if (counter == 1) then
            xval(counter) = 0.0_dp
          else
            ! If we are printing the first point in a path,
            ! It means that the coordinate did not change (otherwise
            ! we would not be printing it). Therefore I do not move
            ! on the x axis, there was a jump in the path here.
            xval(counter) = xval(counter - 1)
          end if
          plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_spts - 1)

          idx_special_points(2*loop_spts - 1) = counter
          xval_special_points(2*loop_spts - 1) = xval(counter)
        end if

        ! This is looping on all points but the first (1 is the first point
        ! after the first in the path)
        do loop_i = 1, kpath_pts(loop_spts)
          counter = counter + 1
          ! Set xval, the x position on the path of the current path
          if (counter == 1) then
            ! This case should never happen but I keep it in for "safety"
            xval(counter) = 0.0_dp
          else
            xval(counter) = xval(counter - 1) + kpath_len(loop_spts)/real(kpath_pts(loop_spts), dp)
          end if
          plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_spts - 1) + &
                                    (kpoint_path%points(:, 2*loop_spts) &
                                     - kpoint_path%points(:, 2*loop_spts - 1))* &
                                    (real(loop_i, dp)/real(kpath_pts(loop_spts), dp))
        end do
        idx_special_points(2*loop_spts) = counter
        xval_special_points(2*loop_spts) = xval(counter)
      end do
      !xval(total_pts)=sum(kpath_len)
      plot_kpoint(:, total_pts) = kpoint_path%points(:, bands_num_spec_points)
    end if
    !
    ! Write out the kpoints in the path
    !
    if (on_root) then
      open (newunit=bndunit, file=trim(seedname)//'_band.kpt', form='formatted')
      write (bndunit, *) total_pts
      do loop_spts = 1, total_pts
        write (bndunit, '(3f12.6,3x,a)') (plot_kpoint(loop_i, loop_spts), loop_i=1, 3), "1.0"
      end do
      close (bndunit)
      !
      ! Write out information on high-symmetry points in the path
      !
      open (newunit=bndunit, file=trim(seedname)//'_band.labelinfo.dat', form='formatted')
      if (kpoint_path%bands_kpt_explicit) then
        do loop_spts = 1, num_spts
          write (bndunit, '(a,3x,I10,3x,4f18.10)') &
            kpoint_path%labels(label_idx_special_points(loop_spts)), &
            idx_special_points(loop_spts), &
            xval_special_points(loop_spts), &
            (plot_kpoint(loop_i, idx_special_points(loop_spts)), loop_i=1, 3)
        end do
      else
        do loop_spts = 1, bands_num_spec_points
          if ((MOD(loop_spts, 2) .eq. 1) .and. &
              (kpath_print_first_point((loop_spts + 1)/2) .eqv. .false.)) cycle
          write (bndunit, '(a,3x,I10,3x,4f18.10)') &
            kpoint_path%labels(loop_spts), &
            idx_special_points(loop_spts), &
            xval_special_points(loop_spts), &
            (plot_kpoint(loop_i, idx_special_points(loop_spts)), loop_i=1, 3)
        end do
      end if
      close (bndunit)
    end if ! on_root
    !
    ! Cut H matrix in real-space
    !
    if (index(band_plot%mode, 'cut') .ne. 0) then
      call plot_cut_hr(real_space_ham, real_lattice, mp_grid, num_wann, &
                       wannier_centres_translated, stdout, error)
      if (allocated(error)) return
    end if
    !
    ! Interpolate the Hamiltonian at each kpoint
    !
    if (ws_region%use_ws_distance) then
      if (index(band_plot%mode, 's-k') .ne. 0) then
        call ws_translate_dist(ws_distance, ws_region, num_wann, &
                               wannier_data%centres, real_lattice, mp_grid, nrpts, &
                               irvec, error, comm, force_recompute=.true.)
        if (allocated(error)) return

      elseif (index(band_plot%mode, 'cut') .ne. 0) then
        call ws_translate_dist(ws_distance, ws_region, num_wann, &
                               wannier_data%centres, real_lattice, mp_grid, nrpts_cut, &
                               irvec_cut, error, comm, force_recompute=.true.)
        if (allocated(error)) return

      else
        call set_error_warn(error, 'Error in plot_interpolate bands: value of bands_plot_mode not recognised', comm)
        return
      end if
    end if
    if (on_root .and. print_output%timing_level > 2) then
      call io_stopwatch_start('plot: interpolate_bands: loop_kpoints', timer)
    end if
    ! [lp] the s-k and cut codes are very similar when use_ws_distance is used, a complete
    !      merge after this point is not impossible
    call comms_array_split(total_pts, counts, displs, comm)
    ! Don't worry about serial run, it is OK to call comms_array_split!
    do loop_kpt = displs(my_node_id) + 1, displs(my_node_id) + counts(my_node_id)
      ! do loop_kpt = 1, total_pts
      ham_kprm = cmplx_0
      !
      if (index(band_plot%mode, 's-k') .ne. 0) then
        do irpt = 1, nrpts
          ! [lp] Shift the WF to have the minimum distance IJ, see also ws_distance.F90
          if (ws_region%use_ws_distance) then
            do j = 1, num_wann
            do i = 1, num_wann
              do ideg = 1, ws_distance%ndeg(i, j, irpt)
                rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), &
                                          real(ws_distance%irdist(:, ideg, i, j, irpt), dp))
                fac = cmplx(cos(rdotk), sin(rdotk), dp) &
                      /real(ndegen(irpt)*ws_distance%ndeg(i, j, irpt), dp)
                ham_kprm(i, j) = ham_kprm(i, j) + fac*ham_r(i, j, irpt)
              end do
            end do
            end do
          else
            ! [lp] Original code, without IJ-dependent shift:
            rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), irvec(:, irpt))
            fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ndegen(irpt), dp)
            ham_kprm = ham_kprm + fac*ham_r(:, :, irpt)
          end if
        end do
        ! end of s-k mode
      elseif (index(band_plot%mode, 'cut') .ne. 0) then
        do irpt = 1, nrpts_cut
          ! [lp] Shift the WF to have the minimum distance IJ, see also ws_distance.F90
          if (ws_region%use_ws_distance) then
            do j = 1, num_wann
            do i = 1, num_wann
              do ideg = 1, ws_distance%ndeg(i, j, irpt)
                rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), &
                                          real(ws_distance%irdist(:, ideg, i, j, irpt), dp))
                fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ws_distance%ndeg(i, j, irpt), dp)
                ham_kprm(i, j) = ham_kprm(i, j) + fac*ham_r_cut(i, j, irpt)
              end do
            end do
            end do
            ! [lp] Original code, without IJ-dependent shift:
          else
            rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), irvec_cut(:, irpt))
            !~[aam] check divide by ndegen?
            fac = cmplx(cos(rdotk), sin(rdotk), dp)
            ham_kprm = ham_kprm + fac*ham_r_cut(:, :, irpt)
          end if ! end of use_ws_distance
        end do
      end if ! end of "cut" mode
      !
      ! Diagonalise H_k (->basis of eigenstates)
      do j = 1, num_wann
        do i = 1, j
          ham_pack(i + ((j - 1)*j)/2) = ham_kprm(i, j)
        end do
      end do
      call ZHPEVX('V', 'A', 'U', num_wann, ham_pack, 0.0_dp, 0.0_dp, 0, 0, -1.0_dp, &
                  nfound, eig_int(1, loop_kpt), U_int, num_wann, cwork, rwork, iwork, ifail, info)
      if (info < 0) then
        write (stdout, '(a,i3,a)') 'THE ', -info, ' ARGUMENT OF ZHPEVX HAD AN ILLEGAL VALUE'
        call set_error_fatal(error, 'Error in plot_interpolate_bands', comm)
        return
      end if
      if (info > 0) then
        write (stdout, '(i3,a)') info, ' EIGENVECTORS FAILED TO CONVERGE'
        call set_error_warn(error, 'Error in plot_interpolate_bands', comm)
        return
      end if
      ! Compute projection onto WF if requested
      if (allocated(band_plot%project)) then
      do loop_w = 1, num_wann
        do loop_p = 1, num_wann
          if (any(band_plot%project == loop_p)) then
            bands_proj(loop_w, loop_kpt) = bands_proj(loop_w, loop_kpt) + &
                                           abs(U_int(loop_p, loop_w))**2
          end if
        end do
      end do
      end if
      !
    end do

    call comms_reduce(eig_int(1, 1), num_wann*total_pts, 'SUM', error, comm)
    call comms_reduce(bands_proj(1, 1), num_wann*total_pts, 'SUM', error, comm)

    if (on_root .and. print_output%timing_level > 2) then
      call io_stopwatch_stop('plot: interpolate_bands: loop_kpoints', timer)
    end if
    !
    ! Interpolation Finished!
    ! Now we write plotting files
    !
    if (on_root) then
      emin = minval(eig_int) - 1.0_dp
      emax = maxval(eig_int) + 1.0_dp

      if (index(band_plot%format, 'gnu') > 0) then
        call plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_points, num_wann)
      end if
      if (index(band_plot%format, 'xmgr') > 0) then
        call plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann, error)
      end if
      write (stdout, '(1x,a,f11.3,a)') &
        'Time to calculate interpolated band structure ', io_time() - time0, ' (sec)'
      write (stdout, *)
    end if

    if (allocated(ham_r_cut)) then
      deallocate (ham_r_cut, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating ham_r_cut in plot_interpolate_bands', comm)
        return
      end if
    end if
    if (allocated(irvec_cut)) then
      deallocate (irvec_cut, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating irvec_cut in plot_interpolate_bands', comm)
        return
      end if
    end if

    if (print_output%timing_level > 1) call io_stopwatch_stop('plot: interpolate_bands', timer)

    if (allocated(idx_special_points)) then
      deallocate (idx_special_points, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating idx_special_points in &
            &plot_interpolate_bands', comm)
        return
      end if
    end if
    if (allocated(xval_special_points)) then
      deallocate (xval_special_points, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating xval_special_points in &
            &plot_interpolate_bands', comm)
        return
      end if
    end if
    if (allocated(label_idx_special_points)) then
      deallocate (label_idx_special_points, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating label_idx_special_points in &
            &plot_interpolate_bands', comm)
        return
      end if
    end if

  contains

    !================================================!
    subroutine plot_cut_hr(real_space_ham, real_lattice, mp_grid, num_wann, &
                           wannier_centres_translated, stdout, error)
      !================================================!
      !
      !!  In real-space picture, ham_r(j,i,k) is an interaction between
      !!  j_th WF at 0 and i_th WF at the lattice point translated
      !!  by matmul(real_lattice(:,:),irvec(:,k))
      !!  We truncate Hamiltonian matrix when
      !!   1) |  r_i(0) - r_j (R) | > dist_cutoff
      !!   2) |  ham_r(i,j,k)     | < hr_cutoff
      !!  while the condition 1) is essential to get a meaningful band structure,
      !!    ( dist_cutoff must be smaller than the shortest distance from
      !!      the center of W-S supercell to the points at the cell boundaries )
      !!  the condition 2) is optional.
      !!
      !!  limitation: when bands_plot_dim .ne. 3
      !!      one_dim_vec must be parallel to one of the cartesian axis
      !!      and perpendicular to the other two primitive lattice vectors
      !================================================!

      use w90_constants, only: dp, cmplx_0, eps8
      use w90_wannier90_types, only: band_plot_type, real_space_ham_type
      use w90_error, only: w90_error_type, set_error_alloc, set_error_warn

      implicit none

      ! arguments
      type(real_space_ham_type), intent(in) :: real_space_ham
      type(w90_error_type), allocatable, intent(out) :: error

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

      integer, intent(in) :: mp_grid(3)
      integer, intent(in) :: num_wann
      integer, intent(in) :: stdout

      ! local variables
      integer :: nrpts_tmp
      integer :: one_dim_vec, two_dim_vec(2)
      integer :: i, j, n1, n2, n3, i1, i2, i3
      real(kind=dp), allocatable :: ham_r_tmp(:, :)
      real(kind=dp), allocatable :: shift_vec(:, :)
      real(kind=dp) :: dist_ij_vec(3)
      real(kind=dp) :: dist_vec(3)
      real(kind=dp) :: dist

      allocate (ham_r_tmp(num_wann, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating ham_r_tmp in plot_cut_hr', comm)
        return
      end if

      irvec_max = maxval(irvec, DIM=2) + 1

      if (real_space_ham%system_dim .ne. 3) then
        ! Find one_dim_vec which is parallel to one_dim_dir
        ! two_dim_vec - the other two lattice vectors
        ! Along the confined directions, take only irvec=0
        j = 0
        do i = 1, 3
          if (abs(abs(real_lattice(real_space_ham%one_dim_dir, i)) &
                  - sqrt(dot_product(real_lattice(:, i), real_lattice(:, i)))) .lt. eps8) then
            one_dim_vec = i
            j = j + 1
          end if
        end do
        if (j .ne. 1) then
          call set_error_warn(error, 'Error: 1-d lattice vector not defined in plot_cut_hr', comm)
          return
        end if
        j = 0
        do i = 1, 3
          if (i .ne. one_dim_vec) then
            j = j + 1
            two_dim_vec(j) = i
          end if
        end do
        if (real_space_ham%system_dim .eq. 1) then
          irvec_max(two_dim_vec(1)) = 0
          irvec_max(two_dim_vec(2)) = 0
        end if
        if (real_space_ham%system_dim .eq. 2) irvec_max(one_dim_vec) = 0
      end if

      nrpts_cut = (2*irvec_max(1) + 1)*(2*irvec_max(2) + 1)*(2*irvec_max(3) + 1)
      allocate (ham_r_cut(num_wann, num_wann, nrpts_cut), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating ham_r_cut in plot_cut_hr', comm)
        return
      end if
      allocate (irvec_cut(3, nrpts_cut), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating irvec_cut in plot_cut_hr', comm)
        return
      end if
      allocate (shift_vec(3, nrpts_cut), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating shift_vec in plot_cut_hr', comm)
        return
      end if

      nrpts_tmp = 0
      do n1 = -irvec_max(1), irvec_max(1)
        do n2 = -irvec_max(2), irvec_max(2)
          loop_n3: do n3 = -irvec_max(3), irvec_max(3)
            do irpt = 1, nrpts
              i1 = mod(n1 - irvec(1, irpt), mp_grid(1))
              i2 = mod(n2 - irvec(2, irpt), mp_grid(2))
              i3 = mod(n3 - irvec(3, irpt), mp_grid(3))
              if (i1 .eq. 0 .and. i2 .eq. 0 .and. i3 .eq. 0) then
                nrpts_tmp = nrpts_tmp + 1
                ham_r_cut(:, :, nrpts_tmp) = ham_r(:, :, irpt)
                irvec_cut(1, nrpts_tmp) = n1
                irvec_cut(2, nrpts_tmp) = n2
                irvec_cut(3, nrpts_tmp) = n3
                cycle loop_n3
              end if
            end do
          end do loop_n3
        end do
      end do

      if (nrpts_tmp .ne. nrpts_cut) then
        write (stdout, '(a)') 'FAILED TO EXPAND ham_r'
        call set_error_warn(error, 'Error in plot_cut_hr', comm)
        return
      end if

      ! AAM: 29/10/2009 Bug fix thanks to Dr Shujun Hu, NIMS, Japan.
      do irpt = 1, nrpts_cut
        ! line below is incorrect for non-orthorhombic cells
        !shift_vec(:,irpt) = matmul(real_lattice(:,:),real(irvec_cut(:,irpt),dp))
        ! line below is the same as calculating
        ! matmul(transpose(real_lattice(:,:)),irvec_cut(:,irpt))
        shift_vec(:, irpt) = matmul(real(irvec_cut(:, irpt), dp), real_lattice(:, :))
      end do

      ! note: dist_cutoff_mode does not necessarily follow bands_plot_dim
      ! e.g. for 1-d system (bands_plot_dim=1) we can still apply 3-d dist_cutoff (dist_cutoff_mode=three_dim)
      if (index(real_space_ham%dist_cutoff_mode, 'one_dim') > 0) then
        do i = 1, num_wann
          do j = 1, num_wann
            dist_ij_vec(real_space_ham%one_dim_dir) = &
              wannier_centres_translated(real_space_ham%one_dim_dir, i) - &
              wannier_centres_translated(real_space_ham%one_dim_dir, j)
            do irpt = 1, nrpts_cut
              dist_vec(real_space_ham%one_dim_dir) = dist_ij_vec(real_space_ham%one_dim_dir) + &
                                                     shift_vec(real_space_ham%one_dim_dir, irpt)
              dist = abs(dist_vec(real_space_ham%one_dim_dir))
              if (dist .gt. real_space_ham%dist_cutoff) &
                ham_r_cut(j, i, irpt) = cmplx_0
            end do
          end do
        end do
      else if (index(real_space_ham%dist_cutoff_mode, 'two_dim') > 0) then
        do i = 1, num_wann
          do j = 1, num_wann
            dist_ij_vec(:) = wannier_centres_translated(:, i) - wannier_centres_translated(:, j)
            do irpt = 1, nrpts_cut
              dist_vec(:) = dist_ij_vec(:) + shift_vec(:, irpt)
              dist_vec(real_space_ham%one_dim_dir) = 0.0_dp
              dist = sqrt(dot_product(dist_vec, dist_vec))
              if (dist .gt. real_space_ham%dist_cutoff) &
                ham_r_cut(j, i, irpt) = cmplx_0
            end do
          end do
        end do
      else
        do i = 1, num_wann
          do j = 1, num_wann
            dist_ij_vec(:) = wannier_centres_translated(:, i) - wannier_centres_translated(:, j)
            do irpt = 1, nrpts_cut
              dist_vec(:) = dist_ij_vec(:) + shift_vec(:, irpt)
              dist = sqrt(dot_product(dist_vec, dist_vec))
              if (dist .gt. real_space_ham%dist_cutoff) &
                ham_r_cut(j, i, irpt) = cmplx_0
            end do
          end do
        end do
      end if

      do irpt = 1, nrpts_cut
        do i = 1, num_wann
          do j = 1, num_wann
            if (abs(ham_r_cut(j, i, irpt)) .lt. real_space_ham%hr_cutoff) &
              ham_r_cut(j, i, irpt) = cmplx_0
          end do
        end do
      end do

      write (stdout, '(/1x,a78)') repeat('-', 78)
      write (stdout, '(1x,4x,a)') &
        'Maximum absolute value of Real-space Hamiltonian at each lattice point'
      write (stdout, '(1x,8x,a62)') repeat('-', 62)
      write (stdout, '(1x,11x,a,11x,a)') 'Lattice point R', 'Max |H_ij(R)|'
      !  output maximum ham_r_cut at each lattice point
      do irpt = 1, nrpts_cut
        ham_r_tmp(:, :) = abs(ham_r_cut(:, :, irpt))
        write (stdout, '(1x,9x,3I5,9x,F12.6)') irvec_cut(:, irpt), maxval(ham_r_tmp)
      end do
      !
      return

    end subroutine plot_cut_hr

    !================================================!
    subroutine plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_points, num_wann)
      !================================================!
      !
      !! Plots the interpolated band structure in gnuplot format
      !
      !================================================!

      use w90_constants, only: dp
      use w90_types, only: kpoint_path_type
      use w90_wannier90_types, only: band_plot_type

      implicit none

      ! arguments
      type(band_plot_type), intent(in) :: band_plot
      type(kpoint_path_type), intent(in) :: kpoint_path
      integer, intent(in) :: num_wann, bands_num_spec_points

      open (newunit=bndunit, file=trim(seedname)//'_band.dat', form='formatted')
      open (newunit=gnuunit, file=trim(seedname)//'_band.gnu', form='formatted')
      !
      ! Gnuplot format
      !
      do i = 1, num_wann
        do nkp = 1, total_pts
          if (allocated(band_plot%project)) then
            write (bndunit, '(3E16.8)') xval(nkp), eig_int(i, nkp), bands_proj(i, nkp)
          else
            write (bndunit, '(2E16.8)') xval(nkp), eig_int(i, nkp)
          end if
        end do
        write (bndunit, *) ' '
      end do
      close (bndunit)
      ! Axis labels
      if (kpoint_path%bands_kpt_explicit) then
        do i = 1, num_spts
          if (i > 1) then
            if (idx_special_points(i) .eq. idx_special_points(i - 1) + 1 .and. &
                kpoint_path%labels(label_idx_special_points(i)) .ne. kpoint_path%labels(label_idx_special_points(i - 1))) then
              ! If two different point indeces are consecutive, the label should be combined
              glabel(i) = TRIM(kpoint_path%labels(label_idx_special_points(i - 1)))//'|'// &
                          TRIM(kpoint_path%labels(label_idx_special_points(i)))
            end if
          else
            glabel(i) = TRIM(kpoint_path%labels(label_idx_special_points(i)))
          end if
        end do
      else
        glabel(1) = TRIM(kpoint_path%labels(1))
        do i = 2, num_paths
          if (kpoint_path%labels(2*(i - 1)) /= kpoint_path%labels(2*(i - 1) + 1)) then
            glabel(i) = TRIM(kpoint_path%labels(2*(i - 1)))//'|'// &
                        TRIM(kpoint_path%labels(2*(i - 1) + 1))
          else
            glabel(i) = TRIM(kpoint_path%labels(2*(i - 1)))
          end if
        end do
        glabel(num_paths + 1) = TRIM(kpoint_path%labels(2*num_paths))
      end if
      ! gnu file
      write (gnuunit, 701) xval(total_pts), emin, emax
      if (kpoint_path%bands_kpt_explicit) then
        do i = 1, num_spts
          write (gnuunit, 705) xval(idx_special_points(i)), emin, xval(idx_special_points(i)), emax
        end do
        write (gnuunit, 702, advance="no") TRIM(glabel(1)), 0.0_dp, &
          (TRIM(glabel(i)), xval(idx_special_points(i)), i=2, num_spts - 1)
        write (gnuunit, 703) TRIM(glabel(num_spts)), xval(idx_special_points(num_spts))
      else
        do i = 1, num_paths - 1
          write (gnuunit, 705) sum(kpath_len(1:i)), emin, sum(kpath_len(1:i)), emax
        end do
        write (gnuunit, 702, advance="no") TRIM(glabel(1)), 0.0_dp, &
          (TRIM(glabel(i + 1)), sum(kpath_len(1:i)), i=1, bands_num_spec_points/2 - 1)
        write (gnuunit, 703) TRIM(glabel(1 + bands_num_spec_points/2)), sum(kpath_len(:))
      end if
      write (gnuunit, *) 'plot ', '"'//trim(seedname)//'_band.dat', '"'
      close (gnuunit)

      if (allocated(band_plot%project)) then
        open (newunit=gnuunit, file=trim(seedname)//'_band_proj.gnu', form='formatted')
        write (gnuunit, '(a)') '#File to plot a colour-mapped Bandstructure'
        write (gnuunit, '(a)') 'set palette defined ( 0 "blue", 3 "green", 6 "yellow", 10 "red" )'
        write (gnuunit, '(a)') 'unset ztics'
        write (gnuunit, '(a)') 'unset key'
        write (gnuunit, '(a)') '# can make pointsize smaller (~0.5). Too small and nothing is &
          &printed'
        write (gnuunit, '(a)') 'set pointsize 0.8'
        write (gnuunit, '(a)') 'set grid xtics'
        write (gnuunit, '(a)') 'set view 0,0'
        write (gnuunit, '(a,f9.5,a)') 'set xrange [0:', xval(total_pts), ']'
        write (gnuunit, '(a,f9.5,a,f9.5,a)') 'set yrange [', emin, ':', emax, ']'
        if (kpoint_path%bands_kpt_explicit) then
          write (gnuunit, 702, advance="no") TRIM(glabel(1)), 0.0_dp, &
            (TRIM(glabel(i)), xval(idx_special_points(i)), i=2, num_spts - 1)
          write (gnuunit, 703) TRIM(glabel(num_spts)), xval(idx_special_points(num_spts))
        else
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(i + 1), sum(kpath_len(1:i)), i=1, bands_num_spec_points/2 - 1)
          write (gnuunit, 703) glabel(1 + bands_num_spec_points/2), sum(kpath_len(:))
        end if

        write (gnuunit, '(a,a,a,a)') 'splot ', '"'//trim(seedname)//'_band.dat', '"', &
          ' u 1:2:3 w p pt 13 palette'
        write (gnuunit, '(a)') '#use the next lines to make a nice figure for a paper'
        write (gnuunit, '(a)') '#set term postscript enhanced eps color lw 0.5 dl 0.5'
        write (gnuunit, '(a)') '#set pointsize 0.275'
      end if
      !
701   format('set style data dots', /, 'set nokey', /, 'set xrange [0:', F8.5, ']', /, &
             'set yrange [', F9.5, ' :', F9.5, ']')
702   format('set xtics (', :20('"', A, '" ', F8.5, ','))
703   format(A, '" ', F8.5, ')')
705   format('set arrow from ', F8.5, ',', F10.5, ' to ', F8.5, ',', F10.5, ' nohead')

    end subroutine plot_interpolate_gnuplot

    !================================================!
    subroutine plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann, error)
      !================================================!
      !
      !! Plots the interpolated band structure in Xmgrace format
      !
      !================================================!

      use w90_io, only: io_date
      use w90_types, only: kpoint_path_type
      use w90_error, only: w90_error_type, set_error_warn

      implicit none

      type(kpoint_path_type), intent(in) :: kpoint_path
      type(w90_error_type), allocatable, intent(out) :: error

      integer, intent(in) :: num_wann, bands_num_spec_points
      character(len=9) :: cdate, ctime

      if (kpoint_path%bands_kpt_explicit) then
        call set_error_warn(error, "bands_kpt_explicit not implemented with xmgrace format", comm)
      end if

      call io_date(cdate, ctime)

      ! Axis labels

      ! Switch any G to Gamma

      do i = 1, bands_num_spec_points
        if (kpoint_path%labels(i) == 'G') then
          ctemp(i) = '\xG\0'
        else
          ctemp(i) = kpoint_path%labels(i)
        end if
      end do

      xlabel(1) = ' '//trim(ctemp(1))//' '
      do i = 2, num_paths
        if (ctemp(2*(i - 1)) /= ctemp(2*(i - 1) + 1)) then
          xlabel(i) = trim(ctemp(2*(i - 1)))//'|'//trim(ctemp(2*(i - 1) + 1))
        else
          xlabel(i) = ctemp(2*(i - 1))
        end if
      end do
      xlabel(num_paths + 1) = ctemp(bands_num_spec_points)

      open (newunit=gnuunit, file=trim(seedname)//'_band.agr', form='formatted')
      !
      ! Xmgrace format
      !
      write (gnuunit, '(a)') '# Grace project file                      '
      write (gnuunit, '(a)') '# written using Wannier90 www.wannier.org '
      write (gnuunit, '(a)') '@version 50113                            '
      write (gnuunit, '(a)') '@page size 792, 612                       '
      write (gnuunit, '(a)') '@page scroll 5%                           '
      write (gnuunit, '(a)') '@page inout 5%                            '
      write (gnuunit, '(a)') '@link page off                            '
      write (gnuunit, '(a)') '@timestamp def "'//cdate//' at '//ctime//'" '
      write (gnuunit, '(a)') '@with g0'
      write (gnuunit, '(a)') '@    world xmin 0.00'
      write (gnuunit, '(a,f10.5)') '@    world xmax ', xval(total_pts)
      write (gnuunit, '(a,f10.5)') '@    world ymin ', emin
      write (gnuunit, '(a,f10.5)') '@    world ymax ', emax
      write (gnuunit, '(a)') '@default linewidth 1.5'
      write (gnuunit, '(a)') '@    xaxis  tick on'
      write (gnuunit, '(a)') '@    xaxis  tick major 1'
      write (gnuunit, '(a)') '@    xaxis  tick major color 1'
      write (gnuunit, '(a)') '@    xaxis  tick major linestyle 3'
      write (gnuunit, '(a)') '@    xaxis  tick major grid on'
      write (gnuunit, '(a)') '@    xaxis  tick spec type both'
      write (gnuunit, '(a,i0)') '@    xaxis  tick spec ', 1 + bands_num_spec_points/2
      write (gnuunit, '(a)') '@    xaxis  tick major 0, 0'
      do i = 1, bands_num_spec_points/2
        write (gnuunit, '(a,i0,a,a)') '@    xaxis  ticklabel ', i - 1, ',', '"'// &
          trim(adjustl(xlabel(i)))//'"'
        write (gnuunit, '(a,i0,a,f10.5)') '@    xaxis  tick major ', i, ' , ', sum(kpath_len(1:i))
      end do
      write (gnuunit, '(a,i0,a)') '@    xaxis  ticklabel ', bands_num_spec_points/2 &
        , ',"'//trim(adjustl(xlabel(1 + bands_num_spec_points/2)))//'"'
      write (gnuunit, '(a)') '@    xaxis  ticklabel char size 1.500000'
      write (gnuunit, '(a)') '@    yaxis  tick major 10'
      write (gnuunit, '(a)') '@    yaxis  label "Band Energy (eV)"'
      write (gnuunit, '(a)') '@    yaxis  label char size 1.500000'
      write (gnuunit, '(a)') '@    yaxis  ticklabel char size 1.500000'
      do i = 1, num_wann
        write (gnuunit, '(a,i0,a)') '@    s', i - 1, ' line color 1'
      end do
      do i = 1, num_wann
        write (gnuunit, '(a,i0)') '@target G0.S', i - 1
        write (gnuunit, '(a)') '@type xy'
        do nkp = 1, total_pts
          write (gnuunit, '(2E16.8)') xval(nkp), eig_int(i, nkp)
        end do
        write (gnuunit, '(a,i0)') '&'
      end do

    end subroutine plot_interpolate_xmgrace

  end subroutine plot_interpolate_bands
plot_interpolate_bands Subroutine

Contents

Source Code

private subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, real_space_ham, ws_region, print_output, recip_lattice, num_wann, wannier_data, ham_r, irvec, ndegen, nrpts, wannier_centres_translated, ws_distance, bands_num_spec_points, stdout, seedname, timer, error, comm)

Uses

Arguments

Calls

Called by

Source Code
