diff --git a/src/plot.F90 b/src/plot.F90 index 10cf05ed..ed866601 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -385,11 +385,12 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, integer :: loop_spts, total_pts, loop_i, nkp, ideg integer :: num_paths, num_spts, ierr integer :: bndunit, gnuunit, loop_w, loop_p - integer :: kpath_pts(bands_num_spec_points/2) + integer, allocatable :: kpath_pts(:) integer, allocatable :: idx_special_points(:) + integer, allocatable :: label_idx_special_points(:) - real(kind=dp) :: kpath_len(bands_num_spec_points/2) 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(:, :) @@ -404,7 +405,7 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, complex(kind=dp), allocatable :: U_int(:, :) complex(kind=dp), allocatable :: cwork(:) - logical :: kpath_print_first_point(bands_num_spec_points/2) + logical, allocatable :: kpath_print_first_point(:) character(len=20), allocatable :: glabel(:) character(len=10), allocatable :: xlabel(:) @@ -474,62 +475,93 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, return endif - 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 - endif - 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 - endif - idx_special_points = -1 - xval_special_points = -1._dp + 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 + endif + 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 + endif + 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 + endif + 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 + endif + 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 + endif + 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 ! - 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 - enddo + 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 - ! 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 + kpath_print_first_point = .false. - 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 + ! 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 + enddo + + ! 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 @@ -553,7 +585,7 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, endif allocate (glabel(num_spts), stat=ierr) if (ierr /= 0) then - call set_error_alloc(error, 'Error in allocating num_spts in plot_interpolate_bands', comm) + call set_error_alloc(error, 'Error in allocating glabel in plot_interpolate_bands', comm) return endif allocate (xlabel(num_spts), stat=ierr) @@ -570,46 +602,90 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, ! ! Find the position of each kpoint in the path ! - 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) - endif - plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_spts - 1) + 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 + endif + 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 + endif + 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 + endif + 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 .and. 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 + exit + end if + end do + end do - idx_special_points(2*loop_spts - 1) = counter - xval_special_points(2*loop_spts - 1) = xval(counter) - end if + 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) + endif + plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_spts - 1) - ! 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) - endif - 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)) + 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) + endif + 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 - 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) + !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 ! @@ -624,15 +700,25 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, ! Write out information on high-symmetry points in the path ! open (newunit=bndunit, file=trim(seedname)//'_band.labelinfo.dat', form='formatted') - 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 + 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) endif ! on_root ! @@ -772,7 +858,7 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, call plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_points, num_wann) endif if (index(band_plot%format, 'xmgr') > 0) then - call plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann) + call plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann, error) endif write (stdout, '(1x,a,f11.3,a)') & 'Time to calculate interpolated band structure ', io_time() - time0, ' (sec)' @@ -812,6 +898,14 @@ subroutine plot_interpolate_bands(mp_grid, real_lattice, band_plot, kpoint_path, return endif endif + 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 + endif + endif contains @@ -1058,24 +1152,46 @@ subroutine plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_point enddo close (bndunit) ! Axis labels - 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)) + if (kpoint_path%bands_kpt_explicit) then + do i = 1, num_spts + if (i > 1 .and. 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))) + 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 - do i = 1, num_paths - 1 - write (gnuunit, 705) sum(kpath_len(1:i)), emin, sum(kpath_len(1:i)), emax - enddo - 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(:)) + 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 + enddo + 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) @@ -1092,9 +1208,15 @@ subroutine plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_point 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, ']' - 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(:)) + 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' @@ -1112,7 +1234,7 @@ subroutine plot_interpolate_gnuplot(band_plot, kpoint_path, bands_num_spec_point end subroutine plot_interpolate_gnuplot !================================================! - subroutine plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann) + subroutine plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann, error) !================================================! ! !! Plots the interpolated band structure in Xmgrace format @@ -1121,13 +1243,20 @@ subroutine plot_interpolate_xmgrace(kpoint_path, bands_num_spec_points, num_wann 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 diff --git a/src/readwrite.F90 b/src/readwrite.F90 index 6b5ecf11..def9374f 100644 --- a/src/readwrite.F90 +++ b/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,59 @@ 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 +1073,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 +4347,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) !================================================! diff --git a/src/types.F90 b/src/types.F90 index f40dec30..02c91d65 100644 --- a/src/types.F90 +++ b/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 diff --git a/src/wannier90_readwrite.F90 b/src/wannier90_readwrite.F90 index 24ad8b61..5da600bf 100644 --- a/src/wannier90_readwrite.F90 +++ b/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,8 @@ 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 +2089,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 +2120,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