User specified list of k-points for the band structure calculation

src/readwrite.F90

@@ -62,6 +62,7 @@ module w90_readwrite
   public :: w90_readwrite_read_kmesh_data
   public :: w90_readwrite_read_kpath
   public :: w90_readwrite_read_kpoints
+  public :: w90_readwrite_read_explicit_kpath
   public :: w90_readwrite_read_lattice
   public :: w90_readwrite_read_mp_grid
   public :: w90_readwrite_read_num_bands
@@ -450,6 +451,56 @@ subroutine w90_readwrite_read_kpath(settings, kpoint_path, ok, bands_plot, error
     endif
   end subroutine w90_readwrite_read_kpath
 
+  subroutine w90_readwrite_read_explicit_kpath(settings, kpoint_path, ok, bands_plot, bohr, error, comm)
+    use w90_error, only: w90_error_type, set_error_input, set_error_alloc
+    implicit none
+    logical, intent(in) :: bands_plot
+    type(kpoint_path_type), intent(inout) :: kpoint_path
+    logical, intent(out) :: ok
+    real(kind=dp), intent(in) :: bohr
+    type(w90_error_type), allocatable, intent(out) :: error
+    type(w90_comm_type), intent(in) :: comm
+    type(settings_type), intent(inout) :: settings
+
+    integer :: i_temp, ierr, bands_num_spec_points
+    logical :: found
+
+    bands_num_spec_points = 0
+    kpoint_path%bands_kpt_explicit = .false.
+    call w90_readwrite_get_block_length(settings, 'explicit_kpath_labels', found, bands_num_spec_points, error, comm)
+    if (allocated(error)) return
+    if (found) then
+      ok = .true.
+      kpoint_path%bands_kpt_explicit = .true.
+!      bands_num_spec_points = i_temp*2
+      if (allocated(kpoint_path%labels)) deallocate (kpoint_path%labels)
+      allocate (kpoint_path%labels(bands_num_spec_points), stat=ierr)
+      if (ierr /= 0) then
+        call set_error_alloc(error, 'Error allocating explicit_kpath labels in w90_wannier90_readwrite_read', comm)
+        return
+      endif
+      if (allocated(kpoint_path%points)) deallocate (kpoint_path%points)
+      allocate (kpoint_path%points(3, bands_num_spec_points), stat=ierr)
+      if (ierr /= 0) then
+        call set_error_alloc(error, 'Error allocating explicit kpoint labels in w90_wannier90_readwrite_read', comm)
+        return
+      endif
+      call w90_readwrite_get_keyword_explicit_kpath(settings, kpoint_path, error, comm)
+      if (allocated(error)) return
+      call w90_readwrite_read_explicit_kpath_points(settings, kpoint_path%bands_kpt_frac, bohr, &
+                                                    error, comm)
+      if (allocated(error)) return
+    else
+      ok = .false.
+    end if
+    ! if (bands_plot) then
+    !   if (kpoint_path%num_points_first_segment < 0) then
+    !     call set_error_input(error, 'Error: bands_num_points must be positive', comm)
+    !     return
+    !   endif
+    ! endif
+  end subroutine w90_readwrite_read_explicit_kpath
+
   subroutine w90_readwrite_read_fermi_energy(settings, found_fermi_energy, fermi_energy_list, &
                                              error, comm)
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
@@ -875,6 +926,58 @@ subroutine w90_readwrite_read_kpoints(settings, pw90_effective_model, kpt_latt,
     endif
   end subroutine w90_readwrite_read_kpoints
 
+  subroutine w90_readwrite_read_explicit_kpath_points(settings, kpt_latt, bohr, &
+                                                      error, comm)
+    use w90_error, only: w90_error_type, set_error_input, set_error_alloc, set_error_dealloc
+    implicit none
+
+    ! arguments
+    real(kind=dp), allocatable, intent(out) :: kpt_latt(:, :)
+    real(kind=dp), intent(in) :: bohr
+    type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
+
+    ! local variables
+    real(kind=dp), allocatable :: kpt_cart(:, :)
+    integer :: ierr, num_kpts
+    logical :: found
+
+    ! pw90_effective_model ignores kpt_cart
+    ! this routine allocates the intent(out) kpt_latt
+
+    call w90_readwrite_get_block_length(settings, 'explicit_kpath', found, num_kpts, error, comm)
+
+    ierr = 0
+
+    allocate (kpt_latt(3, num_kpts), stat=ierr)
+    if (ierr /= 0) then
+      call set_error_alloc(error, 'Error allocating kpt_latt in w90_readwrite_read_explicit_kpath', comm)
+      return
+    endif
+
+    allocate (kpt_cart(3, num_kpts), stat=ierr)
+    if (ierr /= 0) then
+      call set_error_alloc(error, 'Error allocating kpt_cart in w90_readwrite_read_explicit_kpath', comm)
+      return
+    endif
+
+    call w90_readwrite_get_keyword_block(settings, 'explicit_kpath', found, num_kpts, 3, bohr, error, &
+                                         comm, r_value=kpt_cart)
+    if (allocated(error)) return
+    if (.not. found) then
+      call set_error_input(error, 'Error: Found explicit_kpath_labels but there is no explicit_kpath block', comm)
+      return
+    endif
+    kpt_latt = kpt_cart
+
+    deallocate (kpt_cart, stat=ierr)
+    if (ierr /= 0) then
+      call set_error_dealloc(error, 'Error deallocating kpt_cart in w90_readwrite_read_explicit_kpath', comm)
+      return
+    endif
+  end subroutine w90_readwrite_read_explicit_kpath_points
+
   subroutine w90_readwrite_read_lattice(settings, real_lattice, bohr, error, comm)
     use w90_error, only: w90_error_type, set_error_input
     implicit none
@@ -969,6 +1072,8 @@ subroutine w90_readwrite_clear_keywords(settings, comm)
     call clear_block(settings, 'dis_spheres', error, comm)
     call clear_block(settings, 'kpoint_path', error, comm)
     call clear_block(settings, 'kpoints', error, comm)
+    call clear_block(settings, 'explicit_kpath_labels', error, comm)
+    call clear_block(settings, 'explicit_kpath', error, comm)
     call clear_block(settings, 'nnkpts', error, comm)
     call clear_block(settings, 'projections', error, comm)
     call clear_block(settings, 'slwf_centres', error, comm)
@@ -4241,6 +4346,101 @@ subroutine w90_readwrite_get_keyword_kpath(settings, kpoint_path, error, comm)
     return
   end subroutine w90_readwrite_get_keyword_kpath
 
+  subroutine w90_readwrite_get_keyword_explicit_kpath(settings, kpoint_path, error, comm)
+    !================================================!
+    !
+    !!  Fills the explicit_kpath_labels data block
+    !
+    !================================================!
+    use w90_error, only: w90_error_type, set_error_input
+
+    implicit none
+
+    type(kpoint_path_type), intent(inout) :: kpoint_path
+    type(w90_error_type), allocatable, intent(out) :: error
+    type(w90_comm_type), intent(in) :: comm
+    type(settings_type), intent(inout) :: settings
+
+    character(len=20) :: keyword
+    integer           :: ic, in, ins, ine, loop, inner_loop, i, line_e, line_s, counter
+    logical           :: found_e, found_s
+    character(len=maxlen) :: dummy, end_st, start_st
+
+    keyword = "explicit_kpath_labels"
+
+    found_s = .false.
+    found_e = .false.
+
+    start_st = 'begin '//trim(keyword)
+    end_st = 'end '//trim(keyword)
+
+    do loop = 1, settings%num_lines
+      ins = index(settings%in_data(loop), trim(keyword))
+      if (ins == 0) cycle
+      in = index(settings%in_data(loop), 'begin')
+      if (in == 0 .or. in > 1) cycle
+      line_s = loop
+      if (found_s) then
+        call set_error_input(error, 'Error: Found '//trim(start_st)//' more than once in input file', comm)
+        return
+      endif
+      found_s = .true.
+    end do
+
+    do loop = 1, settings%num_lines
+      ine = index(settings%in_data(loop), trim(keyword))
+      if (ine == 0) cycle
+      in = index(settings%in_data(loop), 'end')
+      if (in == 0 .or. in > 1) cycle
+      line_e = loop
+      if (found_e) then
+        call set_error_input(error, 'Error: Found '//trim(end_st)//' more than once in input file', comm)
+        return
+      endif
+      found_e = .true.
+    end do
+
+    if (found_s .and. .not. found_e) then
+      call set_error_input(error, 'Error: Found '//trim(start_st)//' but no '//trim(end_st)//' in input file', comm)
+      return
+    end if
+
+    if (found_s .and. found_e) then
+      if (line_e <= line_s) then
+        call set_error_input(error, 'Error: '//trim(end_st)//' comes before '//trim(start_st)//' in input file', comm)
+        return
+      endif
+    else
+      return !just not found
+    end if
+
+    counter = 0
+    do loop = line_s + 1, line_e - 1
+
+      counter = counter + 1
+      dummy = settings%in_data(loop)
+      read (dummy, *, err=240, end=240) kpoint_path%labels(counter), (kpoint_path%points(i, counter), i=1, 3)
+    end do
+
+    ! Upper case bands labels (eg, x --> X)
+    if (allocated(kpoint_path%labels)) then
+      do loop = 1, size(kpoint_path%labels)
+        do inner_loop = 1, len(kpoint_path%labels(loop))
+          ic = ichar(kpoint_path%labels(loop) (inner_loop:inner_loop))
+          if ((ic .ge. ichar('a')) .and. (ic .le. ichar('z'))) &
+            kpoint_path%labels(loop) (inner_loop:inner_loop) = char(ic + ichar('Z') - ichar('z'))
+        enddo
+      enddo
+    endif
+
+    settings%in_data(line_s:line_e) (1:maxlen) = ' '
+
+    return
+
+240 call set_error_input(error, 'w90_readwrite_get_keyword_kpath: Problem reading explicit kpath '//trim(dummy), comm)
+    return
+  end subroutine w90_readwrite_get_keyword_explicit_kpath
+
   !================================================!
   subroutine clear_block(settings, keyword, error, comm)
     !================================================!

src/types.F90

@@ -237,6 +237,8 @@ module w90_types
     integer :: num_points_first_segment = 100
     character(len=20), allocatable :: labels(:)
     real(kind=dp), allocatable :: points(:, :)
+    logical :: bands_kpt_explicit ! use user provided list of kpoints for bands kpath
+    real(kind=dp), allocatable ::bands_kpt_frac(:, :) ! explicit bands kpoints in fractional coordinate
   end type kpoint_path_type
 
   type settings_data

src/wannier90_readwrite.F90

@@ -224,6 +224,7 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
 
     ! local variables
     logical :: has_kpath
+    logical :: has_explicit_kpath
     integer :: num_exclude_bands
     logical :: found_fermi_energy
     logical :: disentanglement
@@ -292,11 +293,16 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
                                     error, comm)
       if (allocated(error)) return
 
+      if (.not. has_kpath) then
+        call w90_readwrite_read_explicit_kpath(settings, kpoint_path, has_explicit_kpath, w90_calculation%bands_plot, &
+                                               bohr, error, comm)
+        if (allocated(error)) return
+      endif
       call w90_wannier90_readwrite_read_plot_info(settings, wvfn_read, error, comm)
       if (allocated(error)) return
 
       call w90_wannier90_readwrite_read_band_plot(settings, band_plot, num_wann, has_kpath, &
-                                                  w90_calculation%bands_plot, error, comm)
+                                                  has_explicit_kpath, w90_calculation%bands_plot, error, comm)
       if (allocated(error)) return
 
       call w90_wannier90_readwrite_read_wann_plot(settings, wann_plot, num_wann, &
@@ -1075,7 +1081,7 @@ end subroutine w90_wannier90_readwrite_read_plot_info
 
   !================================================!
   subroutine w90_wannier90_readwrite_read_band_plot(settings, band_plot, num_wann, has_kpath, &
-                                                    bands_plot, error, comm)
+                                                    has_explicit_kpath, bands_plot, error, comm)
     !================================================!
     ! Plotting
     !================================================!
@@ -1085,6 +1091,7 @@ subroutine w90_wannier90_readwrite_read_band_plot(settings, band_plot, num_wann,
     integer, intent(in) :: num_wann
     logical, intent(in) :: bands_plot
     logical, intent(in) :: has_kpath
+    logical, intent(in) :: has_explicit_kpath
     type(band_plot_type), intent(inout) :: band_plot
     type(settings_type), intent(inout) :: settings
     type(w90_comm_type), intent(in) :: comm
@@ -1123,8 +1130,9 @@ subroutine w90_wannier90_readwrite_read_band_plot(settings, band_plot, num_wann,
       endif
     endif
 
-    if (.not. has_kpath .and. bands_plot) then
-      call set_error_input(error, 'A bandstructure plot has been requested but there is no kpoint_path block', comm)
+    if ((.not. has_kpath) .and. (.not. has_explicit_kpath) .and. bands_plot) then
+      call set_error_input(error, &
+                           'A bandstructure plot has been requested but there is no kpoint_path or explicit_kpath block', comm)
       return
     endif
 
@@ -2082,10 +2090,16 @@ subroutine w90_wannier90_readwrite_write(atom_data, band_plot, dis_control, dis_
           write (stdout, '(1x,a46,10x,L8,13x,a1)') '|  Plotting interpolated bandstructure       :', w90_calculation%bands_plot, '|'
           bands_num_spec_points = 0
           if (allocated(kpoint_path%labels)) bands_num_spec_points = size(kpoint_path%labels)
-          write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Number of K-path sections                :', &
-            bands_num_spec_points/2, '|'
-          write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Divisions along first K-path section     :', &
-            kpoint_path%num_points_first_segment, '|'
+          if (kpoint_path%bands_kpt_explicit) then
+            write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Number of high-symmetry points           :', bands_num_spec_points, '|'
+            write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Total number of points along K-path      :', &
+              size(kpoint_path%bands_kpt_frac, 2), '|'
+          else
+            write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Number of K-path sections                :', &
+              bands_num_spec_points/2, '|'
+            write (stdout, '(1x,a46,10x,I8,13x,a1)') '|   Divisions along first K-path section     :', &
+              kpoint_path%num_points_first_segment, '|'
+          end if
           write (stdout, '(1x,a46,10x,a8,13x,a1)') '|   Output format                            :', &
             trim(band_plot%format), '|'
           write (stdout, '(1x,a46,10x,a8,13x,a1)') '|   Output mode                              :', &
@@ -2107,15 +2121,26 @@ subroutine w90_wannier90_readwrite_write(atom_data, band_plot, dis_control, dis_
               trim(real_space_ham%dist_cutoff_mode), '|'
           endif
           write (stdout, '(1x,a78)') '*----------------------------------------------------------------------------*'
-          write (stdout, '(1x,a78)') '|   K-space path sections:                                                   |'
+          if (kpoint_path%bands_kpt_explicit) then
+            write (stdout, '(1x,a78)') '|   K-space path high symmetry points:                                       |'
+          else
+            write (stdout, '(1x,a78)') '|   K-space path sections:                                                   |'
+          end if
           if (bands_num_spec_points == 0) then
             write (stdout, '(1x,a78)') '|     None defined                                                           |'
           else
-            do loop = 1, bands_num_spec_points, 2
-              write (stdout, '(1x,a10,1x,a5,1x,3F7.3,5x,a3,1x,a5,1x,3F7.3,3x,a1)') '|    From:', &
-                kpoint_path%labels(loop), (kpoint_path%points(i, loop), i=1, 3), &
-                'To:', kpoint_path%labels(loop + 1), (kpoint_path%points(i, loop + 1), i=1, 3), '|'
-            end do
+            if (kpoint_path%bands_kpt_explicit) then
+              do loop = 1, bands_num_spec_points
+                write (stdout, '(1x,a5,a5,1x,3F7.3,a46)') '|    ', kpoint_path%labels(loop), &
+                  (kpoint_path%points(i, loop), i=1, 3), '                                             |'
+              end do
+            else
+              do loop = 1, bands_num_spec_points, 2
+                write (stdout, '(1x,a10,1x,a5,1x,3F7.3,5x,a3,1x,a5,1x,3F7.3,3x,a1)') '|    From:', &
+                  kpoint_path%labels(loop), (kpoint_path%points(i, loop), i=1, 3), &
+                  'To:', kpoint_path%labels(loop + 1), (kpoint_path%points(i, loop + 1), i=1, 3), '|'
+              end do
+            end if
           end if
           write (stdout, '(1x,a78)') '*----------------------------------------------------------------------------*'
         end if

test-suite/tests/jobconfig

@@ -34,6 +34,27 @@ program = WANNIER90_BANDS_PLOT
 inputs_args = ('silicon.win', '')
 output = silicon_band.dat
 
+# Silicon, 4 valence bands + 4 conduction bands; interpolated bandstructure on the explicitly set kpath,
+# checking the labelinfo.dat file
+[testw90_example03_explicit_kpath_labelinfo]
+program = WANNIER90_LABELINFO
+inputs_args = ('silicon.win', '')
+output = silicon_band.labelinfo.dat
+
+# Silicon, 4 valence bands + 4 conduction bands; interpolated bandstructure on the explicitly set kpath,
+# checking the band.dat file
+[testw90_example03_explicit_kpath_bands_plot]
+program = WANNIER90_BANDS_PLOT
+inputs_args = ('silicon.win', '')
+output = silicon_band.dat
+
+# Silicon, 4 valence bands + 4 conduction bands; interpolated bandstructure on the explicitly set kpath,
+# checking the band.kpt file
+[testw90_example03_explicit_kpath_band_kpt]
+program = WANNIER90_BAND_KPT
+inputs_args = ('silicon.win', '')
+output = silicon_band.kpt
+
 # Copper, states around the Fermi level; Fermi surface
 [testw90_example04]
 program = WANNIER90_WOUT_OK

test-suite/tests/testw90_example03_explicit_kpath_band_kpt/.gitignore

@@ -0,0 +1,7 @@
+silicon_band.dat
+silicon_band.kpt
+silicon_band.labelinfo.dat
+silicon_band.gnu
+silicon.wout
+*.mmn
+*.chk