calc_n_double Function

    function calc_n_double(n_orbs, n_alpha, n_beta) result(n_double)
        ! routine to calculate the number of determinants with different
        ! number of doubly occupied sites..
        integer, intent(in) :: n_orbs, n_alpha, n_beta
        integer, allocatable :: n_double(:)
#ifdef DEBUG_
        character(*), parameter :: this_routine = "calc_n_double"
#endif

        real(dp) :: n_first
        integer :: n_max, n_min, i
        ! the number of possible open shell determinants is:
        ! for the half-filled ms=0 case, it is how often we can distribute
        ! n/2 electrons in n orbitals (since the rest then gets filled up
        ! by the other spin electrons -> nchoosek(n, n/2)
        ! for ms/=0 but still at half-filling it i still the binomial
        ! coefficient of one spin-type

        ! off-half filling, we get more then one possible distribution of
        ! the electrons of the other spin for each distribution of the
        ! first spin type..
        ! so we get some product of binomial of the remaining possible
        ! distributions!

        ! todo: write those formulas up!

        ! todo: a way to create it would be like in my matlab code to
        ! create the heisenberg basis states.. there i even did it in
        ! a ordered fashion.

        ! the number of ways to distribute alpha spins on the lattice
        ! maybe the number of determinants with one double occupancy is
        ! also of interest.. (those are the single excitations from the
        ! fully open shell dets(in the case off-half-filling it could a
        ! single excitation can also lead to another fully-open-shell det, but
        ! these are already covered.
        ! one beta electron MUST reside in a orbital with a alpha spin (thats
        ! just nOccAlpha) and the rest (nOccBeta-1) is on an empty orbital

        ! just for the fun of it: i could calculate the size of each of
        ! those space(1 double occ, 2 double occ.. etc.

        ! only do that for less equal half-filling!

        n_max = max(n_alpha, n_beta)
        n_min = min(n_alpha, n_beta)

        n_first = choose_i64(n_orbs, n_max)

        allocate(n_double(n_min + 1))

        do i = 0, n_min
            n_double(i + 1) = int(n_first * choose_i64(n_max, i) * &
                                  choose_i64(n_orbs - n_max, n_min - i))
        end do

        ! make sure that the sum of basis states is the whole hilber space
        ASSERT(sum(n_double) == choose_i64(n_orbs, n_alpha) * choose_i64(n_orbs, n_beta))

    end function calc_n_double