Main routine
subroutine k_path
!! Main routine
use w90_comms
use w90_constants, only: dp, cmplx_0, cmplx_i, twopi, eps8
use w90_io, only: io_error, io_file_unit, seedname, &
io_time, io_stopwatch, stdout
use w90_utility, only: utility_diagonalize
use w90_postw90_common, only: pw90common_fourier_R_to_k
use w90_parameters, only: num_wann, kpath_task, &
bands_num_spec_points, bands_label, &
kpath_bands_colour, nfermi, fermi_energy_list, &
berry_curv_unit, shc_alpha, shc_beta, shc_gamma, kubo_adpt_smr
use w90_get_oper, only: get_HH_R, HH_R, get_AA_R, get_BB_R, get_CC_R, &
get_FF_R, get_SS_R, get_SHC_R
use w90_spin, only: spin_get_nk
use w90_berry, only: berry_get_imf_klist, berry_get_imfgh_klist, &
berry_get_shc_klist
use w90_constants, only: bohr
integer, dimension(0:num_nodes - 1) :: counts, displs
integer :: i, j, n, num_paths, num_spts, loop_path, loop_kpt, &
total_pts, counter, loop_i, dataunit, gnuunit, pyunit, &
my_num_pts
real(kind=dp) :: ymin, ymax, kpt(3), spn_k(num_wann), &
imf_k_list(3, 3, nfermi), img_k_list(3, 3, nfermi), &
imh_k_list(3, 3, nfermi), Morb_k(3, 3), &
range, zmin, zmax
real(kind=dp) :: shc_k_band(num_wann), shc_k_fermi(nfermi)
real(kind=dp), allocatable, dimension(:) :: 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(:)
! 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(kpath_task, 'bands') > 0
plot_curv = index(kpath_task, 'curv') > 0
plot_morb = index(kpath_task, 'morb') > 0
plot_shc = index(kpath_task, 'shc') > 0
if (on_root) then
if (plot_shc .or. (plot_bands .and. kpath_bands_colour == 'shc')) then
! not allowed to use adpt smr, since adpt smr needs berry_kmesh,
! see line 1837 of berry.F90
if (kubo_adpt_smr) call io_error( &
'Error: Must use fixed smearing when plotting spin Hall conductivity')
end if
if (plot_shc) then
if (nfermi == 0) then
call io_error('Error: must specify Fermi energy')
else if (nfermi /= 1) then
call io_error('Error: kpath plot only accept one Fermi energy, ' &
//'use fermi_energy instead of fermi_energy_min')
end if
end if
end if
call k_path_print_info(plot_bands, plot_curv, plot_morb, plot_shc)
! Set up the needed Wannier matrix elements
call get_HH_R
if (plot_curv .or. plot_morb) call get_AA_R
if (plot_morb) then
call get_BB_R
call get_CC_R
endif
if (plot_shc .or. (plot_bands .and. kpath_bands_colour == 'shc')) then
call get_AA_R
call get_SS_R
call get_SHC_R
endif
if (plot_bands .and. kpath_bands_colour == 'spin') call get_SS_R
if (on_root) then
! 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)
! (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)
! Partition set of k-points into junks
call comms_array_split(total_pts, counts, displs);
!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)
! 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 (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(kpt, HH_R, HH, 0)
call utility_diagonalize(HH, num_wann, my_eig(:, loop_kpt), UU)
!
! Color-code energy bands with the spin projection along the
! chosen spin quantization axis
!
if (kpath_bands_colour == 'spin') then
call spin_get_nk(kpt, spn_k)
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 (kpath_bands_colour == 'shc') then
call berry_get_shc_klist(kpt, shc_k_band=shc_k_band)
my_color(:, loop_kpt) = shc_k_band
end if
end if
if (plot_morb) then
call berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list)
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(kpt, imf_k_list)
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(kpt, shc_k_fermi=shc_k_fermi)
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)
if (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)
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)
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)
end do
end if
if (plot_shc) then
allocate (shc(total_pts))
call comms_gatherv(my_shc, my_num_pts, shc, counts, displs)
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)
!
dataunit = io_file_unit()
open (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. berry_curv_unit == 'bohr2') curv = curv/bohr**2
if (plot_bands .and. kpath_bands_colour == 'shc') then
if (berry_curv_unit == 'bohr2') color = color/bohr**2
end if
if (plot_shc) then
if (berry_curv_unit == 'bohr2') shc = shc/bohr**2
end if
! Axis labels
!
glabel(1) = ' '//bands_label(1)//' '
do i = 2, num_paths
if (bands_label(2*(i - 1)) /= bands_label(2*(i - 1) + 1)) then
glabel(i) = bands_label(2*(i - 1))//'/'//bands_label(2*(i - 1) + 1)
else
glabel(i) = ' '//bands_label(2*(i - 1))//' '
end if
end do
glabel(num_spts) = ' '//bands_label(bands_num_spec_points)//' '
! 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
!
dataunit = io_file_unit()
file_name = trim(seedname)//'-bands.dat'
write (stdout, '(/,3x,a)') file_name
open (dataunit, file=file_name, form='formatted')
do i = 1, num_wann
do loop_kpt = 1, total_pts
if (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
gnuunit = io_file_unit()
file_name = trim(seedname)//'-bands.gnu'
write (stdout, '(/,3x,a)') file_name
open (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 (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 (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 (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
!
pyunit = io_file_unit()
file_name = trim(seedname)//'-bands.py'
write (stdout, '(/,3x,a)') file_name
open (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 (kpath_bands_colour == 'spin' &
.or. kpath_bands_colour == 'shc') write (pyunit, '(a)') "z=data[:,2]"
if (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 (kpath_bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (kpath_bands_colour == 'spin' .or. &
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 (kpath_bands_colour == 'spin' .or. 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
dataunit = io_file_unit()
file_name = trim(seedname)//'-curv.dat'
write (stdout, '(/,3x,a)') file_name
open (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
!
gnuunit = io_file_unit()
file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (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
!
pyunit = io_file_unit()
file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.py'
write (stdout, '(/,3x,a)') file_name
open (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 (berry_curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ $\AA^2$ ]')"
else if (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
dataunit = io_file_unit()
file_name = trim(seedname)//'-morb.dat'
write (stdout, '(/,3x,a)') file_name
open (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
!
gnuunit = io_file_unit()
file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (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
!
pyunit = io_file_unit()
file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.py'
write (stdout, '(/,3x,a)') file_name
open (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
dataunit = io_file_unit()
file_name = trim(seedname)//'-shc.dat'
write (stdout, '(/,3x,a)') file_name
open (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
!
gnuunit = io_file_unit()
file_name = trim(seedname)//'-shc'//'.gnu'
write (stdout, '(/,3x,a)') file_name
open (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
!
pyunit = io_file_unit()
file_name = trim(seedname)//'-shc'//'.py'
write (stdout, '(/,3x,a)') file_name
open (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 (berry_curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + shc_alpha)//achar(119 + shc_beta) &
//"}^{spin"//achar(119 + shc_gamma) &
//"}(\mathbf{k})$ [ $\AA^2$ ]')"
else if (berry_curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + shc_alpha)//achar(119 + shc_beta) &
//"}^{spin"//achar(119 + shc_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
!
pyunit = io_file_unit()
file_name = trim(seedname)//'-bands+shc'//'.py'
write (stdout, '(/,3x,a)') file_name
open (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 (kpath_bands_colour == 'spin' .or. kpath_bands_colour == 'shc') &
write (pyunit, '(a)') "z=data[:,2]"
if (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 (kpath_bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (kpath_bands_colour == 'spin' &
.or. 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 (berry_curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$log_{10}|\Omega_{" &
//achar(119 + shc_alpha)//achar(119 + shc_beta) &
//"}^{spin"//achar(119 + shc_gamma) &
//"}(\mathbf{k})|$ [ $\AA^2$ ]')"
else if (berry_curv_unit == 'bohr2') then
write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
//achar(119 + shc_alpha)//achar(119 + shc_beta) &
//"}^{spin"//achar(119 + shc_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
pyunit = io_file_unit()
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 (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 (kpath_bands_colour == 'spin') write (pyunit, '(a)') "z=data[:,2]"
if (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 (kpath_bands_colour == 'none') then
write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
else if (kpath_bands_colour == 'spin' .or. 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 (berry_curv_unit == 'ang2') then
write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
//"(\mathbf{k})$ [ $\AA^2$ ]')"
else if (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")')
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