pick_random_unocc_impurity Function

    function pick_random_unocc_impurity(nI, ms, pgen) result(i)
        integer, intent(in) :: nI(nel)
        integer, intent(in) :: ms
        real(dp), intent(inout) :: pgen
        integer :: i

        real(dp) :: r
        integer :: pool(nImp)
        integer :: j, nEmpty
        logical :: tEmpty(nImp)

        ! First, get the unoccupied impurity sites
        do j = 1, nImp
            tEmpty(j) = (all(nI /= ImpuritySites(j)) .and. G1(ImpuritySites(j))%Ms == ms)
        end do
        pool = pack(ImpuritySites, tEmpty, vector=ImpuritySites)
        nEmpty = count(tEmpty)

        r = genrand_real2_dSFMT()
        i = pool(int(r * nEmpty) + 1)

        ! Adjust pgen
        pgen = pgen * 1.0 / nEmpty

    end function pick_random_unocc_impurity