postw90 Program

program postw90
  !! The postw90 program
  use w90_constants, only: dp, eps6
  use w90_parameters
  use w90_io

  use w90_kmesh
  use w90_comms, only: on_root, num_nodes, comms_setup, comms_end, comms_bcast, comms_barrier
  use w90_postw90_common, only: pw90common_wanint_setup, pw90common_wanint_get_kpoint_file, &
    pw90common_wanint_param_dist, pw90common_wanint_data_dist

  ! These modules deal with the interpolation of specific physical properties
  !
  use w90_dos
  use w90_berry
  use w90_gyrotropic
  use w90_spin
  use w90_kpath
  use w90_kslice

  use w90_boltzwann
  use w90_geninterp

  implicit none

  integer       :: nkp, len_seedname
  logical       :: have_gamma
  real(kind=dp) :: time0, time1, time2
  character(len=9) :: stat, pos
  logical :: wpout_found, werr_found, dryrun
  character(len=50) :: prog

  ! Put some descriptive comments here
  !
  call comms_setup

  library = .false.
  ispostw90 = .true.

  if (on_root) then
    time0 = io_time()
    prog = 'postw90'
    call io_commandline(prog, dryrun)
    len_seedname = len(seedname)
  end if
  call comms_bcast(len_seedname, 1)
  call comms_bcast(seedname, len_seedname)
  call comms_bcast(dryrun, 1)

  if (on_root) then
    ! If an error file (generated by postw90) exists, I delete it
    ! Note: I do it only for the error file generated from the root node;
    ! If error files generated by other nodes exist, I don't do anything for them
    inquire (file=trim(seedname)//'.node_00000.werr', exist=werr_found)
    if (werr_found) then
      stdout = io_file_unit()
      open (unit=stdout, file=trim(seedname)//'.node_00000.werr', status='old', position='append')
      close (stdout, status='delete')
    end if

    inquire (file=trim(seedname)//'.wpout', exist=wpout_found)
    if (wpout_found) then
      stat = 'old'
    else
      stat = 'replace'
    endif
    pos = 'append'

    stdout = io_file_unit()
    open (unit=stdout, file=trim(seedname)//'.wpout', status=trim(stat), position=trim(pos))

    call param_write_header
    if (num_nodes == 1) then
#ifdef MPI
      write (stdout, '(/,1x,a)') 'Running in serial (with parallel executable)'
#else
      write (stdout, '(/,1x,a)') 'Running in serial (with serial executable)'
#endif
    else
      write (stdout, '(/,1x,a,i3,a/)') &
        'Running in parallel on ', num_nodes, ' CPUs'
    endif
  end if

  ! Read onto the root node all the input parameters from seendame.win,
  ! as well as the energy eigenvalues on the ab-initio q-mesh from seedname.eig
  !
  if (on_root) then
    call param_read
    call param_postw90_write
    time1 = io_time()
    write (stdout, '(1x,a25,f11.3,a)') &
      'Time to read parameters  ', time1 - time0, ' (sec)'

    if (.not. effective_model) then
      ! Check if the q-mesh includes the gamma point
      !
      have_gamma = .false.
      do nkp = 1, num_kpts
        if (all(abs(kpt_latt(:, nkp)) < eps6)) have_gamma = .true.
      end do
      if (.not. have_gamma) write (stdout, '(1x,a)') &
        'Ab-initio does not include Gamma. Interpolation may be incorrect!!!'
      !
      ! Need nntot, wb, and bk to evaluate WF matrix elements of
      ! the position operator in reciprocal space. Also need
      ! nnlist to compute the additional matrix elements entering
      ! the orbital magnetization
      !
      call kmesh_get
      time2 = io_time()
      write (stdout, '(1x,a25,f11.3,a)') &
        'Time to get kmesh        ', time2 - time1, ' (sec)'
    endif

    ! GP, May 10, 2012: for the moment I leave this commented
    ! since we need first to tune that routine so that it doesn't
    ! print the memory information related to wannier90.x.
    ! Note that the code for the memory estimation for the
    ! Boltzwann routine is already there.
    !       call param_memory_estimate
  end if

  if (dryrun) then
    if (on_root) then
      write (stdout, *) ' '
      write (stdout, *) '                       ==============================='
      write (stdout, *) '                                   DRYRUN             '
      write (stdout, *) '                       No problems found with win file'
      write (stdout, *) '                       ==============================='
    endif
    stop
  endif

  ! We now distribute a subset of the parameters to the other nodes
  !
  call pw90common_wanint_param_dist

  if (.not. effective_model) then
    !
    ! Read files seedname.chk (overlap matrices, unitary matrices for
    ! both disentanglement and maximal localization, etc.)
    !
    if (on_root) call param_read_chkpt()
    !
    ! Distribute the information in the um and chk files to the other nodes
    !
    ! Ivo: For interpolation purposes we do not need u_matrix_opt and
    !      u_matrix separately, only their product v_matrix, and this
    !      is what is distributed now
    !
    call pw90common_wanint_data_dist
    !
  end if

  ! Read list of k-points in irreducible BZ and their weights
  !
  ! Should this be done on root node only?
  !
  if (wanint_kpoint_file) call pw90common_wanint_get_kpoint_file

  ! Setup a number of common variables for all interpolation tasks
  !
  call pw90common_wanint_setup

  if (on_root) then
    time1 = io_time()
    write (stdout, '(/1x,a25,f11.3,a)') &
      'Time to read and process .chk    ', time1 - time2, ' (sec)'
  endif
  !
  ! Now perform one or more of the following tasks

  ! ---------------------------------------------------------------
  ! Density of states calculated using a uniform interpolation mesh
  ! ---------------------------------------------------------------
  !
  if (dos .and. index(dos_task, 'dos_plot') > 0) call dos_main

! find_fermi_level commented for the moment in dos.F90
!  if(dos .and. index(dos_task,'find_fermi_energy')>0) call find_fermi_level

  ! --------------------------------------------------------------------
  ! Bands, Berry curvature, or orbital magnetization plot along a k-path
  ! --------------------------------------------------------------------
  !
  if (kpath) call k_path

  ! ---------------------------------------------------------------------------
  ! Bands, Berry curvature, or orbital magnetization plot on a slice in k-space
  ! ---------------------------------------------------------------------------
  !
  if (kslice) call k_slice

  ! --------------------
  ! Spin magnetic moment
  ! --------------------
  !
  if (spin_moment) call spin_get_moment

  ! -------------------------------------------------------------------
  ! dc Anomalous Hall conductivity and eventually (if 'mcd' string also
  ! present in addition to 'ahe', e.g., 'ahe+mcd') dichroic optical
  ! conductivity, both calculated on the same (adaptively-refined) mesh
  ! -------------------------------------------------------------------
  !
  ! ---------------------------------------------------------------
  ! Absorptive dichroic optical conductivity & JDOS on uniform mesh
  ! ---------------------------------------------------------------
  !
  ! -----------------------------------------------------------------
  ! Absorptive ordinary optical conductivity & JDOS on a uniform mesh
  ! -----------------------------------------------------------------
  !
  ! -----------------------------------------------------------------
  ! Orbital magnetization
  ! -----------------------------------------------------------------
  !
  if (berry) call berry_main
  ! -----------------------------------------------------------------
  ! Boltzmann transport coefficients (BoltzWann module)
  ! -----------------------------------------------------------------
  !
  if (on_root) then
    time1 = io_time()
  endif

  if (geninterp) call geninterp_main

  if (boltzwann) call boltzwann_main

  if (gyrotropic) call gyrotropic_main

  if (on_root .and. boltzwann) then
    time2 = io_time()
    write (stdout, '(/1x,a,f11.3,a)') &
      'Time for BoltzWann (Boltzmann transport) ', time2 - time1, ' (sec)'
  endif

  ! I put a barrier here before calling the final time printing routines,
  ! just to be sure that all processes have arrived here.
  call comms_barrier
  if (on_root) then
    write (stdout, '(/,1x,a25,f11.3,a)') &
      'Total Execution Time     ', io_time(), ' (sec)'
    if (timing_level > 0) call io_print_timings()
    write (stdout, *)
    write (stdout, '(/,1x,a)') 'All done: postw90 exiting'
    close (stdout)
  end if

  call comms_end

end program postw90