get_spin_density_neighbors Function

    function get_spin_density_neighbors(ilut, orb) result(spin_dens_neighbors)
        ! function to get the spin-density of the neighboring orbitals
        ! n_{i,\up} - n_{i,\down}
        integer(n_int), intent(in) :: ilut(0:NIfTot)
        integer, intent(in) :: orb
        real(dp) :: spin_dens_neighbors
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_spin_density_neighbors"
#endif
        integer, allocatable :: neighbors(:)
        integer :: i

        ASSERT(associated(lat))

        spin_dens_neighbors = 0.0_dp

        ! orb is given in spatial orbials
        neighbors = lat%get_neighbors(orb)
        do i = 1, size(neighbors)
            ! what do we degine as up spin?? have to be sure here
            ! lets take alpha
            if (IsOcc(ilut, 2 * neighbors(i) - 1)) spin_dens_neighbors = spin_dens_neighbors - 1.0_dp
            if (IsOcc(ilut, 2 * neighbors(i))) spin_dens_neighbors = spin_dens_neighbors + 1.0_dp
        end do

        spin_dens_neighbors = spin_dens_neighbors / 2.0_dp

    end function get_spin_density_neighbors