wann_omega_gamma Subroutine

  subroutine wann_omega_gamma(m_w, csheet, sheet, rave, r2ave, rave2, wann_spread)
    !==================================================================!
    !                                                                  !
    !   Calculate the Wannier Function spread                          !
    !                                                                  !
    !===================================================================
    use w90_parameters, only: num_wann, nntot, wbtot, wb, bk, &
      omega_invariant, timing_level
    use w90_io, only: io_error, io_stopwatch

    implicit none

    real(kind=dp), intent(in)  :: m_w(:, :, :)
    complex(kind=dp), intent(in)  :: csheet(:, :, :)
    real(kind=dp), intent(in)  :: sheet(:, :, :)
    real(kind=dp), intent(out) :: rave(:, :)
    real(kind=dp), intent(out) :: r2ave(:)
    real(kind=dp), intent(out) :: rave2(:)
    type(localisation_vars), intent(out)  :: wann_spread

    !local variables
    real(kind=dp) :: summ, brn
    real(kind=dp), allocatable :: m_w_nn2(:)
    integer :: ind, nn, m, n, iw, rn, cn, ierr

    if (timing_level > 1) call io_stopwatch('wann: omega_gamma', 1)

    allocate (m_w_nn2(num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating m_w_nn2 in wann_omega_gamma')

    if (nntot .eq. 3) then
      do nn = 1, nntot
        rn = 2*nn - 1
        cn = 2*nn
        do n = 1, num_wann
          ln_tmp(n, nn, 1) = atan2(m_w(n, n, cn), m_w(n, n, rn))
        end do
      end do
    else
      do nn = 1, nntot
        rn = 2*nn - 1
        cn = 2*nn
        do n = 1, num_wann
          ln_tmp(n, nn, 1) = aimag(log(csheet(n, nn, 1)*cmplx(m_w(n, n, rn), m_w(n, n, cn), dp))) &
                             - sheet(n, nn, 1)
        end do
      end do
    endif

    rave = 0.0_dp
    do iw = 1, num_wann
      do ind = 1, 3
        do nn = 1, nntot
          rave(ind, iw) = rave(ind, iw) - wb(nn)*bk(ind, nn, 1) &
                          *ln_tmp(iw, nn, 1)
        enddo
      enddo
    enddo

    rave2 = 0.0_dp
    do iw = 1, num_wann
      rave2(iw) = sum(rave(:, iw)*rave(:, iw))
    enddo

    m_w_nn2 = 0.0_dp
    r2ave = wbtot
    do iw = 1, num_wann
      do nn = 1, nntot
        rn = 2*nn - 1
        cn = 2*nn
        m_w_nn2(iw) = m_w_nn2(iw) + m_w(iw, iw, rn)**2 + m_w(iw, iw, cn)**2
        r2ave(iw) = r2ave(iw) + wb(nn)*ln_tmp(iw, nn, 1)**2
      enddo
      r2ave(iw) = r2ave(iw) - m_w_nn2(iw)
    enddo

    if (first_pass) then
      summ = 0.0_dp
      do nn = 1, nntot
        rn = 2*nn - 1
        cn = 2*nn
        do m = 1, num_wann
          do n = 1, num_wann
            summ = summ + m_w(n, m, rn)**2 + m_w(n, m, cn)**2
          enddo
        enddo
      enddo
      wann_spread%om_i = wbtot*real(num_wann, dp) - summ
      first_pass = .false.
    else
      wann_spread%om_i = omega_invariant
    endif

    wann_spread%om_od = wbtot*real(num_wann, dp) - sum(m_w_nn2(:)) - wann_spread%om_i

    if (nntot .eq. 3) then
      wann_spread%om_d = 0.0_dp
    else
      wann_spread%om_d = 0.0_dp
      do nn = 1, nntot
        do n = 1, num_wann
          brn = sum(bk(:, nn, 1)*rave(:, n))
          wann_spread%om_d = wann_spread%om_d + wb(nn) &
                             *(ln_tmp(n, nn, 1) + brn)**2
        enddo
      enddo
    end if

    wann_spread%om_tot = wann_spread%om_i + wann_spread%om_d + wann_spread%om_od

    deallocate (m_w_nn2, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating m_w_nn2 in wann_omega_gamma')

    if (timing_level > 1) call io_stopwatch('wann: omega_gamma', 2)

    return

  end subroutine wann_omega_gamma