plot_fermi_surface Subroutine

private subroutine plot_fermi_surface(fermi_energy_list, recip_lattice, fermi_surface_plot, num_wann, ham_r, irvec, ndegen, nrpts, timing_level, stdout, seedname, timer, error, comm)

Uses

  • proc~~plot_fermi_surface~~UsesGraph proc~plot_fermi_surface plot_fermi_surface module~w90_constants w90_constants proc~plot_fermi_surface->module~w90_constants module~w90_error w90_error proc~plot_fermi_surface->module~w90_error module~w90_io w90_io proc~plot_fermi_surface->module~w90_io module~w90_types w90_types proc~plot_fermi_surface->module~w90_types module~w90_wannier90_types w90_wannier90_types proc~plot_fermi_surface->module~w90_wannier90_types module~w90_comms w90_comms module~w90_error->module~w90_comms module~w90_error_base w90_error_base module~w90_error->module~w90_error_base module~w90_io->module~w90_constants module~w90_types->module~w90_constants module~w90_wannier90_types->module~w90_constants module~w90_comms->module~w90_constants module~w90_comms->module~w90_error_base

Prepares a Xcrysden bxsf file to view the fermi surface

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), allocatable :: fermi_energy_list(:)
real(kind=dp), intent(in) :: recip_lattice(3,3)
type(fermi_surface_plot_type), intent(in) :: fermi_surface_plot
integer, intent(in) :: num_wann
complex(kind=dp), intent(in) :: ham_r(:,:,:)
integer, intent(in) :: irvec(:,:)
integer, intent(in) :: ndegen(:)
integer, intent(in) :: nrpts
integer, intent(in) :: timing_level
integer, intent(in) :: stdout
character(len=50), intent(in) :: seedname
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~plot_fermi_surface~~CallsGraph proc~plot_fermi_surface plot_fermi_surface proc~io_date io_date proc~plot_fermi_surface->proc~io_date proc~io_stopwatch_start io_stopwatch_start proc~plot_fermi_surface->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~plot_fermi_surface->proc~io_stopwatch_stop proc~io_time io_time proc~plot_fermi_surface->proc~io_time proc~set_error_alloc set_error_alloc proc~plot_fermi_surface->proc~set_error_alloc proc~set_error_fatal set_error_fatal proc~plot_fermi_surface->proc~set_error_fatal proc~set_error_warn set_error_warn proc~plot_fermi_surface->proc~set_error_warn zhpevx zhpevx proc~plot_fermi_surface->zhpevx proc~comms_sync_error comms_sync_error proc~set_error_alloc->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_alloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~set_error_warn->proc~comms_sync_error proc~set_error_warn->proc~set_base_error

Called by

proc~~plot_fermi_surface~~CalledByGraph proc~plot_fermi_surface plot_fermi_surface proc~plot_main plot_main proc~plot_main->proc~plot_fermi_surface proc~w90_plot w90_plot proc~w90_plot->proc~plot_main program~wannier wannier program~wannier->proc~w90_plot

Source Code

  subroutine plot_fermi_surface(fermi_energy_list, recip_lattice, fermi_surface_plot, num_wann, &
                                ham_r, irvec, ndegen, nrpts, timing_level, stdout, seedname, &
                                timer, error, comm)
    !================================================!
    !
    !!  Prepares a Xcrysden bxsf file to view the fermi surface
    !
    !================================================!

    use w90_constants, only: dp, cmplx_0, twopi
    use w90_io, only: io_date, io_time, io_stopwatch_start, io_stopwatch_stop
    use w90_wannier90_types, only: fermi_surface_plot_type
    use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal, set_error_warn
    use w90_types, only: timer_list_type

    implicit none

    ! arguments
    type(fermi_surface_plot_type), intent(in)   :: fermi_surface_plot
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm
    complex(kind=dp), intent(in) :: ham_r(:, :, :)
    character(len=50), intent(in)  :: seedname
    real(kind=dp), allocatable, intent(in)      :: fermi_energy_list(:)
    real(kind=dp), intent(in) :: recip_lattice(3, 3)
    integer, intent(in) :: irvec(:, :)
    integer, intent(in) :: ndegen(:)
    integer, intent(in) :: nrpts
    integer, intent(in) :: num_wann
    integer, intent(in) :: stdout
    integer, intent(in) :: timing_level

    ! local variables
    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
    integer              :: fermi_n
    character(len=9)     :: cdate, ctime

    if (timing_level > 1) call io_stopwatch_start('plot: fermi_surface', timer)
    time0 = io_time()
    write (stdout, *)
    write (stdout, '(1x,a)') 'Calculating Fermi surface'
    write (stdout, *)

    fermi_n = 0
    if (allocated(fermi_energy_list)) fermi_n = size(fermi_energy_list)
    if (fermi_n > 1) then
      call set_error_alloc(error, "Error in plot: nfermi>1. Set the fermi level " &
                           //"using the input parameter 'fermi_level'", comm)
      return
    end if

    allocate (ham_pack((num_wann*(num_wann + 1))/2), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ham_pack plot_fermi_surface', comm)
      return
    end if
    allocate (ham_kprm(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ham_kprm plot_fermi_surface', comm)
      return
    end if
    allocate (U_int(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating U_int in plot_fermi_surface', comm)
      return
    end if
    allocate (cwork(2*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cwork in plot_fermi_surface', comm)
      return
    end if
    allocate (rwork(7*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating rwork in plot_fermi_surface', comm)
      return
    end if
    allocate (iwork(5*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating iwork in plot_fermi_surface', comm)
      return
    end if
    allocate (ifail(num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ifail in plot_fermi_surface', comm)
      return
    end if

    npts_plot = (fermi_surface_plot%num_points + 1)**3
    allocate (eig_int(num_wann, npts_plot), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating eig_int in plot_fermi_surface', comm)
      return
    end if
    eig_int = 0.0_dp
    U_int = (0.0_dp, 0.0_dp)

    ikp = 0
    do loop_x = 1, fermi_surface_plot%num_points + 1
      do loop_y = 1, fermi_surface_plot%num_points + 1
        do loop_z = 1, fermi_surface_plot%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_plot%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)
            end do
          end do
          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 set_error_fatal(error, 'Error in plot_fermi_surface', comm)
            return
          end if
          if (info > 0) then
            write (stdout, '(i3,a)') info, ' EIGENVECTORS FAILED TO CONVERGE'
            call set_error_warn(error, 'Error in plot_fermi_surface', comm)
            return
          end if
        end do
      end do
    end do

    call io_date(cdate, ctime)
    open (newunit=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_plot%num_points + 1, fermi_surface_plot%num_points + 1, &
      fermi_surface_plot%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)
      end do
    end do
    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_stop('plot: fermi_surface', timer)

    return

  end subroutine plot_fermi_surface