Prepares a Xcrysden bxsf file to view the fermi surface
subroutine plot_fermi_surface
!===========================================================!
! !
!! Prepares a Xcrysden bxsf file to view the fermi surface
! !
!===========================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i, twopi
use w90_io, only: io_error, stdout, io_file_unit, seedname, &
io_date, io_time, io_stopwatch
use w90_parameters, only: num_wann, fermi_surface_num_points, timing_level, &
recip_lattice, nfermi, fermi_energy_list
use w90_hamiltonian, only: irvec, nrpts, ndegen, ham_r
implicit none
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), allocatable :: eig_int(:, :)
real(kind=dp) :: rdotk, time0
integer, allocatable :: iwork(:), ifail(:)
integer :: loop_x, loop_y, loop_z, INFO, ikp, i, j, ierr
integer :: irpt, nfound, npts_plot, loop_kpt, bxsf_unit
character(len=9) :: cdate, ctime
!
if (timing_level > 1) call io_stopwatch('plot: fermi_surface', 1)
time0 = io_time()
write (stdout, *)
write (stdout, '(1x,a)') 'Calculating Fermi surface'
write (stdout, *)
!
if (nfermi > 1) call io_error("Error in plot: nfermi>1. Set the fermi level " &
//"using the input parameter 'fermi_level'")
!
allocate (ham_pack((num_wann*(num_wann + 1))/2), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ham_pack plot_fermi_surface')
allocate (ham_kprm(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ham_kprm plot_fermi_surface')
allocate (U_int(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating U_int in plot_fermi_surface')
allocate (cwork(2*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cwork in plot_fermi_surface')
allocate (rwork(7*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rwork in plot_fermi_surface')
allocate (iwork(5*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating iwork in plot_fermi_surface')
allocate (ifail(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ifail in plot_fermi_surface')
!
npts_plot = (fermi_surface_num_points + 1)**3
allocate (eig_int(num_wann, npts_plot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating eig_int in plot_fermi_surface')
eig_int = 0.0_dp
U_int = (0.0_dp, 0.0_dp)
!
ikp = 0
do loop_x = 1, fermi_surface_num_points + 1
do loop_y = 1, fermi_surface_num_points + 1
do loop_z = 1, fermi_surface_num_points + 1
ikp = ikp + 1
ham_kprm = cmplx_0
do irpt = 1, nrpts
rdotk = twopi*real((loop_x - 1)*irvec(1, irpt) + &
(loop_y - 1)*irvec(2, irpt) + (loop_z - 1)* &
irvec(3, irpt), dp)/real(fermi_surface_num_points, dp)
fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ndegen(irpt), dp)
ham_kprm = ham_kprm + fac*ham_r(:, :, irpt)
end do
! 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('N', 'A', 'U', num_wann, ham_pack, 0.0_dp, 0.0_dp, 0, 0, -1.0_dp, &
nfound, eig_int(1, ikp), 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_fermi_surface')
endif
if (info > 0) then
write (stdout, '(i3,a)') info, ' EIGENVECTORS FAILED TO CONVERGE'
call io_error('Error in plot_fermi_surface')
endif
end do
end do
end do
call io_date(cdate, ctime)
bxsf_unit = io_file_unit()
open (bxsf_unit, FILE=trim(seedname)//'.bxsf', STATUS='UNKNOWN', FORM='FORMATTED')
write (bxsf_unit, *) ' BEGIN_INFO'
write (bxsf_unit, *) ' #'
write (bxsf_unit, *) ' # this is a Band-XCRYSDEN-Structure-File'
write (bxsf_unit, *) ' # for Fermi Surface Visualisation'
write (bxsf_unit, *) ' #'
write (bxsf_unit, *) ' # Generated by the Wannier90 code http://www.wannier.org'
write (bxsf_unit, *) ' # On ', cdate, ' at ', ctime
write (bxsf_unit, *) ' #'
write (bxsf_unit, *) ' Fermi Energy:', fermi_energy_list(1)
write (bxsf_unit, *) ' END_INFO'
write (bxsf_unit, *)
write (bxsf_unit, *) ' BEGIN_BLOCK_BANDGRID_3D'
write (bxsf_unit, *) 'from_wannier_code'
write (bxsf_unit, *) ' BEGIN_BANDGRID_3D_fermi'
write (bxsf_unit, *) num_wann
write (bxsf_unit, *) fermi_surface_num_points + 1, fermi_surface_num_points + 1, fermi_surface_num_points + 1
write (bxsf_unit, *) '0.0 0.0 0.0'
write (bxsf_unit, *) (recip_lattice(1, i), i=1, 3)
write (bxsf_unit, *) (recip_lattice(2, i), i=1, 3)
write (bxsf_unit, *) (recip_lattice(3, i), i=1, 3)
do i = 1, num_wann
write (bxsf_unit, *) 'BAND: ', i
do loop_kpt = 1, npts_plot
write (bxsf_unit, '(2E16.8)') eig_int(i, loop_kpt)
enddo
enddo
write (bxsf_unit, *) 'END_BANDGRID_3D'
write (bxsf_unit, *) ' END_BLOCK_BANDGRID_3D'
close (bxsf_unit)
write (stdout, '(1x,a,f11.3,a)') 'Time to calculate interpolated Fermi surface ', io_time() - time0, ' (sec)'
write (stdout, *)
!
if (timing_level > 1) call io_stopwatch('plot: fermi_surface', 2)
!
return
end subroutine plot_fermi_surface