calc_number_of_excitations Function

    function calc_number_of_excitations(n_alpha, n_beta, max_excit, n_orbs) &
        result(n_excits)
        ! i might want a routine which calculates the number of all possible
        ! excitations for each excitation level
        integer, intent(in) :: n_alpha, n_beta, max_excit, n_orbs
        integer :: n_excits(max_excit)

        integer :: n_parallel(max_excit, 2), i, j, k

        ! the number of all possible single excitations, per spin is:

        n_parallel = 0
        n_excits = 0

        do i = 1, max_excit
            n_parallel(i, 1) = calc_n_parallel_excitations(n_alpha, n_orbs, i)
            n_parallel(i, 2) = calc_n_parallel_excitations(n_beta, n_orbs, i)
        end do

        ! i can set this up in a general way..

        do i = 1, max_excit

            ! the 2 parallel spin species always are added
            n_excits(i) = sum(n_parallel(i, :))

            ! and i have to zig-zag loop over the other possible combinations
            ! to lead to this level of excitation..
            j = i - 1
            k = 1

            do while (j >= k)

                n_excits(i) = n_excits(i) + n_parallel(j, 1) * n_parallel(k, 2)

                ! and if j /= k do it the other way around too
                if (j /= k) then
                    n_excits(i) = n_excits(i) + n_parallel(j, 2) * n_parallel(k, 1)

                end if
                ! and in/decrease the counters
                j = j - 1
                k = k + 1

                ! in this way i calculate eg.:
                ! n_ij^ab = n_ij_ab(u + d) + n_i^a(u) * n_i^a(d)
                ! and
                ! n_ijk^abc = n_ijk^abc(u + d) + n_ij^ab(u) * n_i^a(d) and spin-flipped..

            end do
        end do

    end function calc_number_of_excitations