!-*- mode: F90 -*-!
!------------------------------------------------------------!
! Copyright (C) 2026 Wannier Developer Group !
! !
! This library is free software; you can redistribute it !
! and/or modify it under the terms of the GNU Lesser General !
! Public License as published by the Free Software !
! Foundation; either version 2.1 of the License, or (at your !
! option) any later version. !
! !
! This library is distributed in the hope that it will be !
! useful,but WITHOUT ANY WARRANTY; without even the implied !
! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR !
! PURPOSE. See the GNU Lesser General Public License for !
! more details. !
! !
! You should have received a copy of the GNU Lesser General !
! Public License along with this library; if not, see !
! <https://www.gnu.org/licenses/>. !
! !
! The webpage of the Wannier90 code is !
! <https://www.wannier.org>. !
! !
! The Wannier90 code is hosted on GitHub !
! <https://github.com/wannier-developers/wannier90> !
!------------------------------------------------------------!
! !
! w90_kpath: evaluate properties along a path in k-space !
! !
!------------------------------------------------------------!
module w90_kpath
!! Calculates quantities along a specified k-path:
!!
!! - Energy bands (eventually colored by the spin)
!! - Energy bands (colored by Berry curvature-like
!! term of spin Hall conductivity)
!!
!! - (Berry curvature)x(-1) summed over occupied bands
!! - Berry curvature-like term of spin Hall conductivity
!! summed over occupied bands
!!
!! - Integrand of orbital magnetization Morb=LCtil+ICtil
use w90_constants, only: dp
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal, &
set_error_input, set_error_fatal, set_error_file
implicit none
private
public :: k_path
contains
!================================================!
! PUBLIC PROCEDURES
!================================================!
subroutine k_path(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, pw90_kpath, kpt_latt, &
pw90_oper_read, pw90_band_deriv_degen, pw90_spin, ws_region, kpoint_path, &
pw90_spin_hall, print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
BB_R, CC_R, HH_R, SH_R, SHR_R, SR_R, SS_R, SAA_R, SBB_R, v_matrix, u_matrix, &
bohr, eigval, real_lattice, scissors_shift, mp_grid, fermi_n, num_wann, &
num_bands, num_kpts, num_valence_bands, effective_model, have_disentangled, &
seedname, stdout, timer, error, comm)
!================================================!
!
!! Main routine
!
!================================================!
use w90_comms, only: w90_comm_type, mpirank, mpisize, comms_array_split, comms_scatterv, &
comms_gatherv, comms_bcast
use w90_constants, only: dp, eps8
use w90_get_oper, only: get_HH_R, get_AA_R_effective, get_AA_R, get_BB_R, get_CC_R, get_SS_R, get_SHC_R
use w90_io, only: io_time
use w90_postw90_common, only: pw90common_fourier_R_to_k
use w90_types, only: kpoint_path_type, print_output_type, wannier_data_type, &
dis_manifold_type, kmesh_info_type, ws_region_type, ws_distance_type, timer_list_type
use w90_postw90_types, only: pw90_berry_mod_type, pw90_spin_hall_type, pw90_kpath_mod_type, &
pw90_spin_mod_type, pw90_band_deriv_degen_type, pw90_oper_read_type, wigner_seitz_type
use w90_berry, only: berry_get_imf_klist, berry_get_imfgh_klist, berry_get_shc_klist
use w90_spin, only: spin_get_nk
use w90_utility, only: utility_diagonalize, utility_recip_lattice_base
implicit none
! arguments
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
type(kmesh_info_type), intent(in) :: kmesh_info
type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
type(pw90_spin_mod_type), intent(in) :: pw90_spin
type(print_output_type), intent(in) :: print_output
type(ws_region_type), intent(in) :: ws_region
type(kpoint_path_type), intent(in) :: kpoint_path
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
type(w90_comm_type), intent(in) :: comm
type(wannier_data_type), intent(in) :: wannier_data
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(ws_distance_type), intent(inout) :: ws_distance
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :) ! <0|H(r-R)|R>
complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :) ! <0|r_alpha.H(r-R)_beta|R>
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :) !<0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :) !<0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), intent(in) :: v_matrix(:, :, :), u_matrix(:, :, :)
real(kind=dp), intent(in) :: bohr
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
real(kind=dp), intent(in) :: kpt_latt(:, :)
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann, num_bands, num_kpts, num_valence_bands, fermi_n
integer, intent(in) :: stdout
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
! local variables
real(kind=dp) :: recip_lattice(3, 3), volume
integer :: i, j, n, num_paths, num_spts, loop_kpt, &
total_pts, loop_i, dataunit, gnuunit, pyunit, &
my_num_pts
real(kind=dp) :: ymin, ymax, kpt(3), spn_k(num_wann), &
imf_k_list(3, 3, fermi_n), img_k_list(3, 3, fermi_n), &
imh_k_list(3, 3, fermi_n), Morb_k(3, 3), &
range, zmin, zmax
real(kind=dp) :: shc_k_band(num_wann), shc_k_fermi(fermi_n)
real(kind=dp), allocatable :: kpath_len(:)
logical :: plot_bands, plot_curv, plot_morb, plot_shc
character(len=120) :: file_name
complex(kind=dp), allocatable :: HH(:, :)
complex(kind=dp), allocatable :: UU(:, :)
real(kind=dp), allocatable :: xval(:), eig(:, :), my_eig(:, :), &
curv(:, :), my_curv(:, :), &
morb(:, :), my_morb(:, :), &
color(:, :), my_color(:, :), &
plot_kpoint(:, :), my_plot_kpoint(:, :), &
shc(:), my_shc(:)
character(len=3), allocatable :: glabel(:)
integer, allocatable :: counts(:), displs(:)
logical :: on_root = .false.
integer :: my_node_id, num_nodes
my_node_id = mpirank(comm)
num_nodes = mpisize(comm)
allocate (counts(0:num_nodes - 1))
allocate (displs(0:num_nodes - 1))
if (my_node_id == 0) on_root = .true.
! Everything is done on the root node (not worthwhile parallelizing)
! However, we still have to read and distribute the data if we
! are in parallel. So calls to get_oper are done on all nodes at the moment
plot_bands = index(pw90_kpath%task, 'bands') > 0
plot_curv = index(pw90_kpath%task, 'curv') > 0
plot_morb = index(pw90_kpath%task, 'morb') > 0
plot_shc = index(pw90_kpath%task, 'shc') > 0
if (on_root) then
if (plot_shc .or. (plot_bands .and. pw90_kpath%bands_colour == 'shc')) then
! not allowed to use adpt smr, since adpt smr needs berry_kmesh,
! see line 1837 of berry.F90
if (pw90_berry%kubo_smearing%use_adaptive) then
call set_error_input(error, 'Error: Must use fixed smearing when plotting spin Hall conductivity', comm)
return
end if
end if
if (plot_shc) then
if (fermi_n == 0) then
call set_error_input(error, 'Error: must specify Fermi energy', comm)
return
else if (fermi_n /= 1) then
call set_error_input(error, 'Error: kpath plot only accept one Fermi energy, ' &
//'use fermi_energy instead of fermi_energy_min', comm)
return
end if
end if
end if
call k_path_print_info(plot_bands, plot_curv, plot_morb, plot_shc, fermi_energy_list, &
pw90_kpath, pw90_berry%curv_unit, stdout, error, comm)
if (allocated(error)) return
! Set up the needed Wannier matrix elements
call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
stdout, timer, error, comm)
if (allocated(error)) return
if (plot_curv .or. plot_morb) then
if (effective_model) then
call get_AA_R_effective(print_output, AA_R, HH_R, wigner_seitz%nrpts, num_wann, seedname, &
stdout, timer, error, comm)
else
call get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wannier_data, AA_R, &
v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
num_wann, have_disentangled, seedname, stdout, timer, error, comm)
end if
if (allocated(error)) return
end if
if (plot_morb) then
call get_BB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, HH_R, BB_R, v_matrix, &
eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
num_wann, have_disentangled, seedname, stdout, timer, error, comm)
if (allocated(error)) return
call get_CC_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
HH_R, BB_R, CC_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, &
ws_region, num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, &
timer, error, comm)
if (allocated(error)) return
end if
if (plot_shc .or. (plot_bands .and. pw90_kpath%bands_colour == 'shc')) then
if (effective_model) then
call get_AA_R_effective(print_output, AA_R, HH_R, wigner_seitz%nrpts, num_wann, seedname, &
stdout, timer, error, comm)
else
call get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wannier_data, AA_R, &
v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
num_wann, have_disentangled, seedname, stdout, timer, error, comm)
end if
if (allocated(error)) return
call get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, eigval, &
wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
have_disentangled, seedname, stdout, timer, error, comm)
if (allocated(error)) return
call get_SHC_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
pw90_spin_hall, SH_R, SHR_R, SR_R, v_matrix, eigval, scissors_shift, &
wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)
if (allocated(error)) return
end if
if (plot_bands .and. pw90_kpath%bands_colour == 'spin') then
call get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, eigval, &
wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
have_disentangled, seedname, stdout, timer, error, comm)
if (allocated(error)) return
end if
num_paths = 0
if (on_root) then
call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
! Determine the number of k-points (total_pts) as well as
! their reciprocal-lattice coordinates long the path (plot_kpoint)
! and their associated horizontal coordinate for the plot (xval)
call k_path_get_points(num_paths, kpath_len, total_pts, xval, plot_kpoint, kpoint_path, &
recip_lattice, pw90_kpath)
! (paths)
num_spts = num_paths + 1 ! number of path endpoints (special pts)
else
! Dummy allocation for making scatterv work
allocate (plot_kpoint(1, 1))
end if
! Broadcast number of k-points on the path
call comms_bcast(total_pts, 1, error, comm)
if (allocated(error)) return
! Partition set of k-points into junks
! call comms_array_split(total_pts, counts, displs)
call comms_array_split(total_pts, counts, displs, comm)
!kpt_lo = displs(my_node_id)+1
!kpt_hi = displs(my_node_id)+counts(my_node_id)
my_num_pts = counts(my_node_id)
! Distribute coordinates
allocate (my_plot_kpoint(3, my_num_pts))
call comms_scatterv(my_plot_kpoint, 3*my_num_pts, &
plot_kpoint, 3*counts, 3*displs, error, comm)
if (allocated(error)) return
! Value of the vertical coordinate in the actual plots: energy bands
!
if (plot_bands) then
allocate (HH(num_wann, num_wann))
allocate (UU(num_wann, num_wann))
allocate (my_eig(num_wann, my_num_pts))
if (pw90_kpath%bands_colour /= 'none') allocate (my_color(num_wann, my_num_pts))
end if
! Value of the vertical coordinate in the actual plots
!
if (plot_curv) allocate (my_curv(my_num_pts, 3))
if (plot_morb) allocate (my_morb(my_num_pts, 3))
if (plot_shc) allocate (my_shc(my_num_pts))
! Loop over local junk of k-points on the path and evaluate the requested quantities
!
do loop_kpt = 1, my_num_pts
kpt(:) = my_plot_kpoint(:, loop_kpt)
if (plot_bands) then
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, &
HH_R, kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
if (allocated(error)) return
call utility_diagonalize(HH, num_wann, my_eig(:, loop_kpt), UU, error, comm)
if (allocated(error)) return
!
! Color-code energy bands with the spin projection along the
! chosen spin quantization axis
!
if (pw90_kpath%bands_colour == 'spin') then
call spin_get_nk(ws_region, pw90_spin, wannier_data, ws_distance, wigner_seitz, HH_R, &
SS_R, kpt, real_lattice, spn_k, mp_grid, num_wann, error, comm)
if (allocated(error)) return
my_color(:, loop_kpt) = spn_k(:)
!
! The following is needed to prevent bands from disappearing
! when the magnitude of the Wannier interpolated spn_k (very
! slightly) exceeds 1.0 (e.g. in bcc Fe along N--G--H)
!
do n = 1, num_wann
if (my_color(n, loop_kpt) > 1.0_dp - eps8) then
my_color(n, loop_kpt) = 1.0_dp - eps8
else if (my_color(n, loop_kpt) < -1.0_dp + eps8) then
my_color(n, loop_kpt) = -1.0_dp + eps8
end if
end do
else if (pw90_kpath%bands_colour == 'shc') then
call berry_get_shc_klist(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
pw90_band_deriv_degen, ws_region, pw90_spin_hall, print_output, &
wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, SH_R, &
SHR_R, SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, &
kpt, real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm, &
shc_k_band=shc_k_band)
if (allocated(error)) return
my_color(:, loop_kpt) = shc_k_band
end if
end if
if (plot_morb) then
call berry_get_imfgh_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, &
real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm, &
imf_k_list, img_k_list, imh_k_list)
if (allocated(error)) return
Morb_k = img_k_list(:, :, 1) + imh_k_list(:, :, 1) &
- 2.0_dp*fermi_energy_list(1)*imf_k_list(:, :, 1)
Morb_k = -Morb_k/2.0_dp ! differs by -1/2 from Eq.97 LVTS12
my_morb(loop_kpt, 1) = sum(Morb_k(:, 1))
my_morb(loop_kpt, 2) = sum(Morb_k(:, 2))
my_morb(loop_kpt, 3) = sum(Morb_k(:, 3))
end if
if (plot_curv) then
if (.not. plot_morb) then
call berry_get_imf_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, &
real_lattice, imf_k_list, scissors_shift, mp_grid, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm)
if (allocated(error)) return
end if
my_curv(loop_kpt, 1) = sum(imf_k_list(:, 1, 1))
my_curv(loop_kpt, 2) = sum(imf_k_list(:, 2, 1))
my_curv(loop_kpt, 3) = sum(imf_k_list(:, 3, 1))
end if
if (plot_shc) then
call berry_get_shc_klist(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
pw90_band_deriv_degen, ws_region, pw90_spin_hall, print_output, &
wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, SH_R, SHR_R, &
SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, kpt, &
real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm, &
shc_k_fermi=shc_k_fermi)
if (allocated(error)) return
my_shc(loop_kpt) = shc_k_fermi(1)
end if
end do !loop_kpt
! Send results to root process
if (plot_bands) then
allocate (eig(num_wann, total_pts))
call comms_gatherv(my_eig, num_wann*my_num_pts, &
eig, num_wann*counts, num_wann*displs, error, comm)
if (allocated(error)) return
if (pw90_kpath%bands_colour /= 'none') then
allocate (color(num_wann, total_pts))
call comms_gatherv(my_color, num_wann*my_num_pts, &
color, num_wann*counts, num_wann*displs, error, comm)
if (allocated(error)) return
end if
end if
if (plot_curv) then
allocate (curv(total_pts, 3))
do i = 1, 3
call comms_gatherv(my_curv(:, i), my_num_pts, &
curv(:, i), counts, displs, error, comm)
if (allocated(error)) return
end do
end if
if (plot_morb) then
allocate (morb(total_pts, 3))
do i = 1, 3
call comms_gatherv(my_morb(:, i), my_num_pts, &
morb(:, i), counts, displs, error, comm)
if (allocated(error)) return
end do
end if
if (plot_shc) then
allocate (shc(total_pts))
call comms_gatherv(my_shc, my_num_pts, shc, counts, displs, error, comm)
if (allocated(error)) return
end if
if (on_root) then
num_spts = num_paths + 1
allocate (glabel(num_spts))
if (plot_bands) then
!
! Write out the kpoints in the path in a format that can be inserted
! directly in the pwscf input file (the '1.0_dp' in the second column is
! a k-point weight, expected by pwscf)
!
open (newunit=dataunit, file=trim(seedname)//'-path.kpt', form='formatted')
write (dataunit, *) total_pts
do loop_kpt = 1, total_pts
write (dataunit, '(3f12.6,3x,f4.1)') &
(plot_kpoint(loop_i, loop_kpt), loop_i=1, 3), 1.0_dp
end do
close (dataunit)
end if
if (plot_curv .and. pw90_berry%curv_unit == 'bohr2') curv = curv/bohr**2
if (plot_bands .and. pw90_kpath%bands_colour == 'shc') then
if (pw90_berry%curv_unit == 'bohr2') color = color/bohr**2
end if
if (plot_shc) then
if (pw90_berry%curv_unit == 'bohr2') shc = shc/bohr**2
end if
! Axis labels
!
glabel(1) = ' '//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) = kpoint_path%labels(2*(i - 1))//'/'//kpoint_path%labels(2*(i - 1) + 1)
else
glabel(i) = ' '//kpoint_path%labels(2*(i - 1))//' '
end if
end do
glabel(num_spts) = ' '//kpoint_path%labels(num_paths*2)//' '
! Now write the plotting files
write (stdout, '(/,1x,a)') 'Output files:'
if (plot_bands) then
file_name = trim(seedname)//'-path.kpt'
write (stdout, '(/,3x,a)') file_name
!
! Data file
!
file_name = trim(seedname)//'-bands.dat'
write (stdout, '(/,3x,a)') file_name
open (newunit=dataunit, file=file_name, form='formatted')
do i = 1, num_wann
do loop_kpt = 1, total_pts
if (pw90_kpath%bands_colour == 'none') then
write (dataunit, '(2E16.8)') xval(loop_kpt), eig(i, loop_kpt)
else
write (dataunit, '(3E16.8)') xval(loop_kpt), &
eig(i, loop_kpt), color(i, loop_kpt)
end if
end do
write (dataunit, *) ' '
end do
close (dataunit)
end if
if (plot_bands .and. .not. plot_curv .and. .not. plot_morb &
.and. .not. plot_shc) then
!
! Gnuplot script
!
ymin = minval(eig) - 1.0_dp
ymax = maxval(eig) + 1.0_dp
file_name = trim(seedname)//'-bands.gnu'
write (stdout, '(/,3x,a)') file_name
open (newunit=gnuunit, file=file_name, form='formatted')
do i = 1, num_paths - 1
write (gnuunit, 705) sum(kpath_len(1:i)), ymin, &
sum(kpath_len(1:i)), ymax
end do
if (pw90_kpath%bands_colour == 'none') then
write (gnuunit, 701) xval(total_pts), ymin, ymax
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) 'plot ', '"'//trim(seedname)//'-bands.dat', '"'
else if (pw90_kpath%bands_colour == 'spin') then
!
! Only works with gnuplot v4.2 and higher
!
write (gnuunit, 706) xval(total_pts), ymin, ymax
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) &
'set palette defined (-1 "blue", 0 "green", 1 "red")'
write (gnuunit, *) 'set pm3d map'
write (gnuunit, *) 'set zrange [-1:1]'
write (gnuunit, *) 'splot ', '"'//trim(seedname)//'-bands.dat', &
'" with dots palette'
else if (pw90_kpath%bands_colour == 'shc') then
!
! Only works with gnuplot v4.2 and higher
!
write (gnuunit, 706) xval(total_pts), ymin, ymax
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) &
'sgnlog10(x) = abs(x) > 10.0 ? log10(abs(x))*sgn(x) : x/10.0'
! merge: Fortran ternary operator: similar to ? : in C
zmin = minval(color)
zmin = merge(sign(log10(abs(zmin)), zmin), zmin/10.0_dp, abs(zmin) > 10.0_dp)
zmax = maxval(color)
zmax = merge(sign(log10(abs(zmax)), zmax), zmax/10.0_dp, abs(zmax) > 10.0_dp)
write (gnuunit, *) &
'set palette defined (', zmin, ' "blue", 0 "green", ', zmax, ' "red")'
write (gnuunit, *) 'set pm3d map'
write (gnuunit, *) 'set zrange [', zmin, ':', zmax, ']'
write (gnuunit, *) 'splot ', '"'//trim(seedname)//'-bands.dat', &
'" u 1:2:(sgnlog10($3)) with dots palette'
end if
close (gnuunit)
!
! python script
!
file_name = trim(seedname)//'-bands.py'
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-bands.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,1]"
if (pw90_kpath%bands_colour == 'spin' &
.or. pw90_kpath%bands_colour == 'shc') write (pyunit, '(a)') "z=data[:,2]"
if (pw90_kpath%bands_colour == 'shc') write (pyunit, '(a)') &
"z=np.array([np.log10(abs(elem))*np.sign(elem) " &
//"if abs(elem)>10 else elem/10.0 for elem in z])"
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
if (pw90_kpath%bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (pw90_kpath%bands_colour == 'spin' .or. &
pw90_kpath%bands_colour == 'shc') then
write (pyunit, '(a)') &
"pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
end if
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
write (pyunit, '(a)') "pl.ylabel('Energy [eV]')"
if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') then
write (pyunit, '(a)') &
"pl.axes().set_aspect(aspect=0.65*max(x)/(max(y)-min(y)))"
write (pyunit, '(a)') "pl.colorbar(shrink=0.7)"
end if
write (pyunit, '(a)') "outfile = '"//trim(seedname)//"-bands.pdf'"
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end if ! plot_bands .and. .not.plot_curv .and. .not.plot_morb .and. .not. plot_shc
if (plot_curv) then
! It is conventional to plot the negative curvature
curv = -curv
file_name = trim(seedname)//'-curv.dat'
write (stdout, '(/,3x,a)') file_name
open (newunit=dataunit, file=file_name, form='formatted')
do loop_kpt = 1, total_pts
write (dataunit, '(4E16.8)') xval(loop_kpt), &
curv(loop_kpt, :)
end do
write (dataunit, *) ' '
close (dataunit)
end if
if (plot_curv .and. .not. plot_bands) then
do i = 1, 3
!
! gnuplot script
!
file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (newunit=gnuunit, file=file_name, form='formatted')
ymin = minval(curv(:, i))
ymax = maxval(curv(:, i))
range = ymax - ymin
ymin = ymin - 0.02_dp*range
ymax = ymax + 0.02_dp*range
write (gnuunit, 707) xval(total_pts), ymin, ymax
do j = 1, num_paths - 1
write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
sum(kpath_len(1:j)), ymax
end do
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) &
'plot ', '"'//trim(seedname)//'-curv.dat', '" u 1:'//achar(49 + i)
close (gnuunit)
!
! python script
!
file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.py'
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-curv.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') "tick_labels.append('" &
//trim(glabel(j))//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
write (pyunit, '(a)') "pl.plot(x,y,color='k')"
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
if (pw90_berry%curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ $\AA^2$ ]')"
else if (pw90_berry%curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ bohr$^2$ ]')"
end if
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-curv_"//achar(119 + i)//".pdf'"
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end do
end if ! plot_curv .and. .not.plot_bands
if (plot_morb) then
file_name = trim(seedname)//'-morb.dat'
write (stdout, '(/,3x,a)') file_name
open (newunit=dataunit, file=file_name, form='formatted')
do loop_kpt = 1, total_pts
write (dataunit, '(4E16.8)') xval(loop_kpt), morb(loop_kpt, :)
end do
write (dataunit, *) ' '
close (dataunit)
end if
if (plot_morb .and. .not. plot_bands) then
do i = 1, 3
!
! gnuplot script
!
file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (newunit=gnuunit, file=file_name, form='formatted')
ymin = minval(morb(:, i))
ymax = maxval(morb(:, i))
range = ymax - ymin
ymin = ymin - 0.02_dp*range
ymax = ymax + 0.02_dp*range
write (gnuunit, 707) xval(total_pts), ymin, ymax
do j = 1, num_paths - 1
write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
sum(kpath_len(1:j)), ymax
end do
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) &
'plot ', '"'//trim(seedname)//'-morb.dat', '" u 1:'//achar(49 + i)
close (gnuunit)
!
! python script
!
file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.py'
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-morb.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') &
"tick_labels.append('"//trim(glabel(j))//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
write (pyunit, '(a)') "pl.plot(x,y,color='k')"
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
write (pyunit, '(a)') "pl.ylabel(r'$M^{\rm{orb}}_z(\mathbf{k})$" &
//" [ Ry$\cdot\AA^2$ ]')"
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-morb_"//achar(119 + i)//".pdf'"
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end do
end if ! plot_morb .and. .not.plot_bands
if (plot_shc) then
file_name = trim(seedname)//'-shc.dat'
write (stdout, '(/,3x,a)') file_name
open (newunit=dataunit, file=file_name, form='formatted')
do loop_kpt = 1, total_pts
write (dataunit, '(2E16.8)') xval(loop_kpt), &
shc(loop_kpt)
end do
write (dataunit, *) ' '
close (dataunit)
end if
if (plot_shc .and. .not. plot_bands) then
!
! gnuplot script
!
file_name = trim(seedname)//'-shc'//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (newunit=gnuunit, file=file_name, form='formatted')
ymin = minval(shc(:))
ymax = maxval(shc(:))
range = ymax - ymin
ymin = ymin - 0.02_dp*range
ymax = ymax + 0.02_dp*range
write (gnuunit, 707) xval(total_pts), ymin, ymax
do j = 1, num_paths - 1
write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
sum(kpath_len(1:j)), ymax
end do
write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
(glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
write (gnuunit, *) &
'plot ', '"'//trim(seedname)//'-shc.dat', '" u 1:2'
close (gnuunit)
!
! python script
!
file_name = trim(seedname)//'-shc'//'.py'
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') "# uncomment these two lines if you are " &
//"running in non-GUI environment"
write (pyunit, '(a)') "#import matplotlib"
write (pyunit, '(a)') "#matplotlib.use('Agg')"
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-shc.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,1]"
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') "tick_labels.append('" &
//trim(glabel(j))//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
write (pyunit, '(a)') "pl.plot(x,y,color='k')"
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
if (pw90_berry%curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
//"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
//"}(\mathbf{k})$ [ $\AA^2$ ]')"
else if (pw90_berry%curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
//"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
//"}(\mathbf{k})$ [ bohr$^2$ ]')"
end if
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-shc"//".pdf'"
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end if ! plot_shc .and. .not.plot_bands
if (plot_bands .and. plot_shc) then
!
! python script
!
file_name = trim(seedname)//'-bands+shc'//'.py'
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') "# uncomment these two lines if you are " &
//"running in non-GUI environment"
write (pyunit, '(a)') "#import matplotlib"
write (pyunit, '(a)') "#matplotlib.use('Agg')"
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') 'from matplotlib.gridspec import GridSpec'
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
write (pyunit, '(a)') "fig = pl.figure()"
write (pyunit, '(a)') "gs = GridSpec(2,52,hspace=0.00,wspace=1)"
!
! upper panel (energy bands)
!
write (pyunit, '(a)') "axes1 = pl.subplot(gs[0, :-2])"
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-bands.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a,F12.6)') "y=data[:,1]-", fermi_energy_list(1)
if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') &
write (pyunit, '(a)') "z=data[:,2]"
if (pw90_kpath%bands_colour == 'shc') &
write (pyunit, '(a)') "z=np.array([np.log10(abs(elem))*np.sign(elem) " &
//"if abs(elem)>10 else elem/10.0 for elem in z])"
if (pw90_kpath%bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (pw90_kpath%bands_colour == 'spin' &
.or. pw90_kpath%bands_colour == 'shc') then
write (pyunit, '(a)') &
"pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
end if
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([-0.65,0.65]) # Adjust this range as needed"
write (pyunit, '(a)') "pl.plot([tick_locs[0],tick_locs[-1]],[0,0]," &
//"color='black',linestyle='--',linewidth=0.5)"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
write (pyunit, '(a)') "pl.ylabel('Energy$-$E$_F$ [eV]')"
write (pyunit, '(a)') "pl.tick_params(axis='x'," &
//"which='both',bottom='off',top='off',labelbottom='off')"
!
! color bar for upper panel
!
write (pyunit, '(a)') "# Now adding the colorbar"
write (pyunit, '(a)') "cbaxes = pl.subplot(gs[:, -2:])"
write (pyunit, '(a)') "cb = pl.colorbar(cax = cbaxes," &
//"orientation='vertical') "
write (pyunit, '(a)') "#cblim = int(min(abs(max(z)),abs(min(z))))"
write (pyunit, '(a)') "#cb.set_ticks([-cblim,0,cblim])"
write (pyunit, '(a)') "#cblim = [min(z),0,max(z)]"
write (pyunit, '(a)') "#cb.set_ticks([min(z),0,max(z)])"
write (pyunit, '(a)') "#cb.set_ticklabels(['-','0','+'])"
!
! lower panel (SHC)
!
write (pyunit, '(a)') "axes2 = pl.subplot(gs[1, :-2])"
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-shc.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,1]"
write (pyunit, '(a)') &
"y=np.array([np.log10(abs(elem))*np.sign(elem) &
&if abs(elem)>10 else elem/10.0 for elem in y])"
write (pyunit, '(a)') "pl.plot(x,y,color='k')"
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.plot([0,max(x)],[0,0],color='black'," &
//"linestyle='--',linewidth=0.5)"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
if (pw90_berry%curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$log_{10}|\Omega_{" &
//achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
//"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
//"}(\mathbf{k})|$ [ $\AA^2$ ]')"
else if (pw90_berry%curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
//"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
//"}(\mathbf{k})$ [ bohr$^2$ ]')"
end if
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-bands+shc"//".pdf'"
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end if ! plot_bands .and. plot_shc
if (plot_bands .and. (plot_curv .or. plot_morb)) then
!
! python script
!
do i = 1, 3
if (plot_curv) then
file_name = trim(seedname)//'-bands+curv_'//achar(119 + i)//'.py'
else if (plot_morb) then
file_name = trim(seedname)//'-bands+morb_'//achar(119 + i)//'.py'
end if
write (stdout, '(/,3x,a)') file_name
open (newunit=pyunit, file=file_name, form='formatted')
write (pyunit, '(a)') 'import pylab as pl'
write (pyunit, '(a)') 'import numpy as np'
write (pyunit, '(a)') 'from matplotlib.gridspec import GridSpec'
write (pyunit, '(a)') "tick_labels=[]"
write (pyunit, '(a)') "tick_locs=[]"
do j = 1, num_spts
if (trim(glabel(j)) == ' G') then
write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
else
write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
//"'.strip())"
end if
if (j == 1) then
write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
else
write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
sum(kpath_len(1:j - 1)), ")"
end if
end do
write (pyunit, '(a)') "fig = pl.figure()"
write (pyunit, '(a)') "gs = GridSpec(2, 1,hspace=0.00)"
!
! upper panel (energy bands)
!
write (pyunit, '(a)') "axes1 = pl.subplot(gs[0, 0:])"
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-bands.dat')"
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a,F12.6)') "y=data[:,1]-", fermi_energy_list(1)
if (pw90_kpath%bands_colour == 'spin') write (pyunit, '(a)') "z=data[:,2]"
if (pw90_kpath%bands_colour == 'shc') then
write (pyunit, '(a)') "z=data[:,2]"
write (pyunit, '(a)') "z=np.array([np.log10(abs(elem))*np.sign(elem) " &
//"if abs(elem)>10 else elem/10.0 for elem in z])"
end if
if (pw90_kpath%bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') then
write (pyunit, '(a)') &
"pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
end if
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([-0.65,0.65]) # Adjust this range as needed"
write (pyunit, '(a)') "pl.plot([tick_locs[0],tick_locs[-1]],[0,0]," &
//"color='black',linestyle='--',linewidth=0.5)"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
write (pyunit, '(a)') "pl.ylabel('Energy$-$E$_F$ [eV]')"
write (pyunit, '(a)') "pl.tick_params(axis='x'," &
//"which='both',bottom='off',top='off',labelbottom='off')"
!
! lower panel (curvature or orbital magnetization)
!
write (pyunit, '(a)') "axes2 = pl.subplot(gs[1, 0:])"
if (plot_curv) then
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-curv.dat')"
else if (plot_morb) then
write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
"-morb.dat')"
end if
write (pyunit, '(a)') "x=data[:,0]"
write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
write (pyunit, '(a)') "pl.plot(x,y,color='k')"
write (pyunit, '(a)') "pl.xlim([0,max(x)])"
write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
//"max(y)+0.025*(max(y)-min(y))])"
write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
write (pyunit, '(a)') " pl.plot([tick_locs[n],tick_locs[n]]," &
//"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
//"linestyle='-',linewidth=0.5)"
if (plot_curv) then
if (pw90_berry%curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ $\AA^2$ ]')"
else if (pw90_berry%curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ bohr$^2$ ]')"
end if
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-bands+curv_"//achar(119 + i)//".pdf'"
else if (plot_morb) then
write (pyunit, '(a)') "pl.ylabel(r'$M^{\rm{orb}}_z(\mathbf{k})$" &
//" [ Ry$\cdot\AA^2$ ]')"
write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
"-morb_"//achar(119 + i)//".pdf'"
end if
write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
write (pyunit, '(a)') "pl.show()"
end do
end if ! plot_bands .and. plot_curv
end if ! on_root
701 format('set style data dots', /, 'unset key', /, &
'set xrange [0:', F8.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')
702 format('set xtics (', :20('"', A3, '" ', F8.5, ','))
703 format(A3, '" ', F8.5, ')')
!704 format('set palette defined (', F8.5, ' "red", 0 "green", ', F8.5, ' "blue")') !not used
705 format('set arrow from ', F16.8, ',', F16.8, ' to ', F16.8, ',', F16.8, ' nohead')
706 format('unset key', /, &
'set xrange [0:', F9.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')
707 format('set style data lines', /, 'set nokey', /, &
'set xrange [0:', F8.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')
end subroutine k_path
!================================================!
! PRIVATE PROCEDURES !
!================================================!
subroutine k_path_print_info(plot_bands, plot_curv, plot_morb, plot_shc, fermi_energy_list, &
pw90_kpath, berry_curv_unit, stdout, error, comm)
!================================================!
use w90_postw90_types, only: pw90_kpath_mod_type
use w90_comms, only: w90_comm_type, mpirank
! arguments
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: stdout
character(len=*), intent(in) :: berry_curv_unit
logical, intent(in) :: plot_bands, plot_curv, plot_morb, plot_shc
! local variables
integer :: fermi_n
logical :: on_root = .false.
if (mpirank(comm) == 0) on_root = .true.
if (on_root) then
write (stdout, '(/,/,1x,a)') &
'Properties calculated in module k p a t h'
write (stdout, '(1x,a)') &
'------------------------------------------'
if (plot_bands) then
select case (pw90_kpath%bands_colour)
case ("none")
write (stdout, '(/,3x,a)') '* Energy bands in eV'
case ("spin")
write (stdout, '(/,3x,a)') '* Energy bands in eV, coloured by spin'
case ("shc")
write (stdout, '(/,3x,a)') '* Energy bands in eV, coloured by SHC'
end select
end if
fermi_n = 0
if (allocated(fermi_energy_list)) fermi_n = size(fermi_energy_list)
if (plot_curv) then
if (berry_curv_unit == 'ang2') then
write (stdout, '(/,3x,a)') '* Negative Berry curvature in Ang^2'
else if (berry_curv_unit == 'bohr2') then
write (stdout, '(/,3x,a)') '* Negative Berry curvature in Bohr^2'
end if
if (fermi_n /= 1) then
call set_error_input(error, 'Must specify one Fermi level when kpath_task=curv', comm)
return
end if
end if
if (plot_morb) then
write (stdout, '(/,3x,a)') &
'* Orbital magnetization k-space integrand in eV.Ang^2'
if (fermi_n /= 1) then
call set_error_input(error, 'Must specify one Fermi level when kpath_task=morb', comm)
return
end if
end if
if (plot_shc) then
if (berry_curv_unit == 'ang2') then
write (stdout, '(/,3x,a)') '* Berry curvature-like term for' &
//' spin Hall conductivity in Ang^2'
else if (berry_curv_unit == 'bohr2') then
write (stdout, '(/,3x,a)') '* Berry curvature-like term for' &
//' spin Hall conductivity in Bohr^2'
end if
if (fermi_n /= 1) then
call set_error_input(error, 'Must specify one Fermi level when kpath_task=shc', comm)
return
end if
end if
end if ! on_root
end subroutine
!================================================!
subroutine k_path_get_points(num_paths, kpath_len, total_pts, xval, plot_kpoint, kpoint_path, &
recip_lattice, pw90_kpath)
!================================================!
! Determine the number of k-points (total_pts) as well as !
! their reciprocal-lattice coordinates long the path (plot_kpoint) !
! and their associated horizontal coordinate for the plot (xval) !
!================================================!
use w90_postw90_types, only: pw90_kpath_mod_type
use w90_types, only: kpoint_path_type
use w90_utility, only: utility_metric
! arguments
type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
type(kpoint_path_type), intent(in) :: kpoint_path
integer, intent(out) :: num_paths, total_pts
real(kind=dp), allocatable, intent(out) :: kpath_len(:), xval(:)
real(kind=dp), allocatable, intent(out) :: plot_kpoint(:, :)
real(kind=dp), intent(in) :: recip_lattice(3, 3)
! local variables
integer :: counter, loop_path, loop_i
integer, allocatable :: kpath_pts(:)
real(kind=dp) :: vec(3)
real(kind=dp) :: recip_metric(3, 3)
call utility_metric(recip_lattice, recip_metric)
! Work out how many points there are in the total path, and the
! positions of the special points
!
if (allocated(kpoint_path%labels)) then
!num_paths = kpoint_path%bands_num_spec_points/2 ! number of straight line segments
num_paths = size(kpoint_path%labels)/2
else
num_paths = 0
end if
allocate (kpath_pts(num_paths))
allocate (kpath_len(num_paths))
do loop_path = 1, num_paths
vec = kpoint_path%points(:, 2*loop_path) &
- kpoint_path%points(:, 2*loop_path - 1)
kpath_len(loop_path) = &
sqrt(dot_product(vec, (matmul(recip_metric, vec))))
!
! kpath_pts(loop_path) is the number of points in path number
! loop_path (all segments have the same density of points)
!
if (loop_path == 1) then
kpath_pts(loop_path) = pw90_kpath%num_points
else
kpath_pts(loop_path) = nint(real(pw90_kpath%num_points, dp) &
*kpath_len(loop_path)/kpath_len(1))
end if
end do
total_pts = sum(kpath_pts) + 1
! Reciprocal-lattice coordinates of the k-points along the path
!
allocate (plot_kpoint(3, total_pts))
! Value of the horizontal coordinate in the actual plots (units of
! distance in k-space)
!
allocate (xval(total_pts))
! Find the position of each kpoint along the path
!
counter = 0
do loop_path = 1, num_paths
do loop_i = 1, kpath_pts(loop_path)
counter = counter + 1
if (counter == 1) then
xval(counter) = 0.0_dp
else
xval(counter) = xval(counter - 1) &
+ kpath_len(loop_path)/real(kpath_pts(loop_path), dp)
end if
plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_path - 1) &
+ (kpoint_path%points(:, 2*loop_path) &
- kpoint_path%points(:, 2*loop_path - 1) &
) &
*(real(loop_i - 1, dp)/real(kpath_pts(loop_path), dp))
end do
end do
!
! Last point
!
xval(total_pts) = sum(kpath_len)
plot_kpoint(:, total_pts) = kpoint_path%points(:, num_paths*2)
end subroutine
end module w90_kpath
