Plots the interpolated band structure
subroutine plot_interpolate_bands
!============================================!
! !
!! Plots the interpolated band structure
! !
!============================================!
use w90_constants, only: dp, cmplx_0, cmplx_i, twopi
use w90_io, only: io_error, stdout, io_file_unit, seedname, &
io_time, io_stopwatch
use w90_parameters, only: num_wann, bands_num_points, recip_metric, &
bands_num_spec_points, timing_level, &
bands_spec_points, bands_label, bands_plot_format, &
bands_plot_mode, num_bands_project, bands_plot_project, &
use_ws_distance
use w90_hamiltonian, only: irvec, nrpts, ndegen, ham_r
use w90_ws_distance, only: irdist_ws, wdist_ndeg, &
ws_translate_dist
implicit none
complex(kind=dp), allocatable :: ham_r_cut(:, :, :)
complex(kind=dp), allocatable :: ham_pack(:)
complex(kind=dp) :: fac
complex(kind=dp), allocatable :: ham_kprm(:, :)
complex(kind=dp), allocatable :: U_int(:, :)
complex(kind=dp), allocatable :: cwork(:)
real(kind=dp), allocatable :: rwork(:)
real(kind=dp) :: kpath_len(bands_num_spec_points/2)
integer :: kpath_pts(bands_num_spec_points/2)
logical :: kpath_print_first_point(bands_num_spec_points/2)
real(kind=dp), allocatable :: xval(:)
real(kind=dp), allocatable :: eig_int(:, :), plot_kpoint(:, :)
real(kind=dp), allocatable :: bands_proj(:, :)
real(kind=dp) :: rdotk, vec(3), emin, emax, time0
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
character(len=20), allocatable :: glabel(:)
character(len=10), allocatable :: xlabel(:)
character(len=10), allocatable :: ctemp(:)
integer, allocatable :: idx_special_points(:)
real(kind=dp), allocatable :: xval_special_points(:)
!
if (timing_level > 1) call io_stopwatch('plot: interpolate_bands', 1)
!
time0 = io_time()
write (stdout, *)
write (stdout, '(1x,a)') 'Calculating interpolated band-structure'
write (stdout, *)
!
allocate (ham_pack((num_wann*(num_wann + 1))/2), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ham_pack in plot_interpolate_bands')
allocate (ham_kprm(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ham_kprm in plot_interpolate_bands')
allocate (U_int(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating U_int in plot_interpolate_bands')
allocate (cwork(2*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cwork in plot_interpolate_bands')
allocate (rwork(7*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rwork in plot_interpolate_bands')
allocate (iwork(5*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating iwork in plot_interpolate_bands')
allocate (ifail(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ifail in plot_interpolate_bands')
allocate (idx_special_points(bands_num_spec_points), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating idx_special_points in plot_interpolate_bands')
allocate (xval_special_points(bands_num_spec_points), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating xval_special_points in plot_interpolate_bands')
idx_special_points = -1
xval_special_points = -1._dp
!
! 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((bands_spec_points(:, (i - 1)*2) - bands_spec_points(:, (i - 1)*2 + 1))**2) > 1.e-6) .or. &
(TRIM(bands_label((i - 1)*2)) .ne. TRIM(bands_label((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 = bands_spec_points(:, 2*loop_spts) - bands_spec_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) = bands_num_points
else
kpath_pts(loop_spts) = nint(real(bands_num_points, 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
allocate (plot_kpoint(3, total_pts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating plot_kpoint in plot_interpolate_bands')
allocate (xval(total_pts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating xval in plot_interpolate_bands')
allocate (eig_int(num_wann, total_pts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating eig_int in plot_interpolate_bands')
allocate (bands_proj(num_wann, total_pts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating bands_proj in plot_interpolate_bands')
allocate (glabel(num_spts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating num_spts in plot_interpolate_bands')
allocate (xlabel(num_spts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating xlabel in plot_interpolate_bands')
allocate (ctemp(bands_num_spec_points), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ctemp in plot_interpolate_bands')
eig_int = 0.0_dp; bands_proj = 0.0_dp
!
! 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) = bands_spec_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)
endif
plot_kpoint(:, counter) = bands_spec_points(:, 2*loop_spts - 1) + &
(bands_spec_points(:, 2*loop_spts) - bands_spec_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) = bands_spec_points(:, bands_num_spec_points)
!
! Write out the kpoints in the path
!
bndunit = io_file_unit()
open (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
!
bndunit = io_file_unit()
open (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)') &
bands_label(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
close (bndunit)
!
! Cut H matrix in real-space
!
if (index(bands_plot_mode, 'cut') .ne. 0) call plot_cut_hr()
!
! Interpolate the Hamiltonian at each kpoint
!
if (use_ws_distance) then
if (index(bands_plot_mode, 's-k') .ne. 0) then
call ws_translate_dist(nrpts, irvec, force_recompute=.true.)
elseif (index(bands_plot_mode, 'cut') .ne. 0) then
call ws_translate_dist(nrpts_cut, irvec_cut, force_recompute=.true.)
else
call io_error('Error in plot_interpolate bands: value of bands_plot_mode not recognised')
endif
endif
! [lp] the s-k and cut codes are very similar when use_ws_distance is used, a complete
! mercge after this point is not impossible
do loop_kpt = 1, total_pts
ham_kprm = cmplx_0
!
if (index(bands_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 (use_ws_distance) then
do j = 1, num_wann
do i = 1, num_wann
do ideg = 1, wdist_ndeg(i, j, irpt)
rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), &
real(irdist_ws(:, ideg, i, j, irpt), dp))
fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ndegen(irpt)*wdist_ndeg(i, j, irpt), dp)
ham_kprm(i, j) = ham_kprm(i, j) + fac*ham_r(i, j, irpt)
enddo
enddo
enddo
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)
endif
end do
! end of s-k mode
elseif (index(bands_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 (use_ws_distance) then
do j = 1, num_wann
do i = 1, num_wann
do ideg = 1, wdist_ndeg(i, j, irpt)
rdotk = twopi*dot_product(plot_kpoint(:, loop_kpt), &
real(irdist_ws(:, ideg, i, j, irpt), dp))
fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(wdist_ndeg(i, j, irpt), dp)
ham_kprm(i, j) = ham_kprm(i, j) + fac*ham_r_cut(i, j, irpt)
enddo
enddo
enddo
! [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)
endif ! end of use_ws_distance
end do
endif ! 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)
enddo
enddo
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 io_error('Error in plot_interpolate_bands')
endif
if (info > 0) then
write (stdout, '(i3,a)') info, ' EIGENVECTORS FAILED TO CONVERGE'
call io_error('Error in plot_interpolate_bands')
endif
! Compute projection onto WF if requested
if (num_bands_project > 0) then
do loop_w = 1, num_wann
do loop_p = 1, num_wann
if (any(bands_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
!
! Interpolation Finished!
! Now we write plotting files
!
emin = minval(eig_int) - 1.0_dp
emax = maxval(eig_int) + 1.0_dp
if (index(bands_plot_format, 'gnu') > 0) call plot_interpolate_gnuplot
if (index(bands_plot_format, 'xmgr') > 0) call plot_interpolate_xmgrace
write (stdout, '(1x,a,f11.3,a)') &
'Time to calculate interpolated band structure ', io_time() - time0, ' (sec)'
write (stdout, *)
if (allocated(ham_r_cut)) deallocate (ham_r_cut, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating ham_r_cut in plot_interpolate_bands')
if (allocated(irvec_cut)) deallocate (irvec_cut, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating irvec_cut in plot_interpolate_bands')
!
if (timing_level > 1) call io_stopwatch('plot: interpolate_bands', 2)
!
if (allocated(idx_special_points)) deallocate (idx_special_points, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating idx_special_points in plot_interpolate_bands')
if (allocated(xval_special_points)) deallocate (xval_special_points, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating xval_special_points in plot_interpolate_bands')
contains
!============================================!
subroutine plot_cut_hr()
!============================================!
!
!! 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_io, only: io_error, stdout
use w90_parameters, only: num_wann, mp_grid, real_lattice, &
one_dim_dir, bands_plot_dim, &
hr_cutoff, dist_cutoff, dist_cutoff_mode
use w90_hamiltonian, only: wannier_centres_translated
implicit none
!
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) call io_error('Error in allocating ham_r_tmp in plot_cut_hr')
irvec_max = maxval(irvec, DIM=2) + 1
if (bands_plot_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(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) call io_error('Error: 1-d lattice vector not defined in plot_cut_hr')
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 (bands_plot_dim .eq. 1) then
irvec_max(two_dim_vec(1)) = 0
irvec_max(two_dim_vec(2)) = 0
end if
if (bands_plot_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) call io_error('Error in allocating ham_r_cut in plot_cut_hr')
allocate (irvec_cut(3, nrpts_cut), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating irvec_cut in plot_cut_hr')
allocate (shift_vec(3, nrpts_cut), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating shift_vec in plot_cut_hr')
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 io_error('Error in plot_cut_hr')
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(dist_cutoff_mode, 'one_dim') > 0) then
do i = 1, num_wann
do j = 1, num_wann
dist_ij_vec(one_dim_dir) = &
wannier_centres_translated(one_dim_dir, i) - wannier_centres_translated(one_dim_dir, j)
do irpt = 1, nrpts_cut
dist_vec(one_dim_dir) = dist_ij_vec(one_dim_dir) + shift_vec(one_dim_dir, irpt)
dist = abs(dist_vec(one_dim_dir))
if (dist .gt. dist_cutoff) &
ham_r_cut(j, i, irpt) = cmplx_0
end do
end do
end do
else if (index(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(one_dim_dir) = 0.0_dp
dist = sqrt(dot_product(dist_vec, dist_vec))
if (dist .gt. 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. 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. 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
!============================================!
! !
!! Plots the interpolated band structure in gnuplot format
!============================================!
use w90_constants, only: dp
use w90_io, only: io_file_unit, seedname
use w90_parameters, only: num_wann, bands_num_spec_points, &
bands_label, num_bands_project
implicit none
!
bndunit = io_file_unit()
open (bndunit, file=trim(seedname)//'_band.dat', form='formatted')
gnuunit = io_file_unit()
open (gnuunit, file=trim(seedname)//'_band.gnu', form='formatted')
!
! Gnuplot format
!
do i = 1, num_wann
do nkp = 1, total_pts
if (num_bands_project > 0) 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
enddo
write (bndunit, *) ' '
enddo
close (bndunit)
! Axis labels
glabel(1) = TRIM(bands_label(1))
do i = 2, num_paths
if (bands_label(2*(i - 1)) /= bands_label(2*(i - 1) + 1)) then
glabel(i) = TRIM(bands_label(2*(i - 1)))//'|'//TRIM(bands_label(2*(i - 1) + 1))
else
glabel(i) = TRIM(bands_label(2*(i - 1)))
end if
end do
glabel(num_paths + 1) = TRIM(bands_label(2*num_paths))
! 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(:))
write (gnuunit, *) 'plot ', '"'//trim(seedname)//'_band.dat', '"'
close (gnuunit)
if (num_bands_project > 0) then
gnuunit = io_file_unit()
open (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, ']'
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(:))
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
!============================================!
! !
!! Plots the interpolated band structure in Xmgrace format
!============================================!
use w90_io, only: io_file_unit, seedname, io_date
use w90_parameters, only: num_wann, bands_num_spec_points
implicit none
character(len=9) :: cdate, ctime
call io_date(cdate, ctime)
! Axis labels
! Switch any G to Gamma
do i = 1, bands_num_spec_points
if (bands_label(i) == 'G') then
ctemp(i) = '\xG\0'
else
ctemp(i) = bands_label(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)
gnuunit = io_file_unit()
open (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