calcNonOverlapDouble Subroutine

    subroutine calcNonOverlapDouble(ilut, csf_i, excitInfo, excitations, nExcits, &
                                    posSwitches, negSwitches)
        ! specific subroutine to calculate the non overlap double excitations
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        integer(n_int), intent(out), allocatable :: excitations(:, :)
        integer, intent(out) :: nExcits
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)

        integer(n_int), allocatable :: tempExcits(:, :), tempExcits2(:, :), tmp_excitations(:, :)
        integer :: tmpNum, iEx, tmpNum2, nOpen1, nOpen2, ierr, iEx2, nMax, i, j
        type(ExcitationInformation_t) :: tmpInfo

        ! plan is to first calculate all lower excitation, and for each state
        ! there calculate the top excitation
        ! have to check that matrix element does not get reset!

        ! modifiy excitInfo so it specifies a single excitation would be better,
        ! or just be lazy for now and call the non excitInfo one.. todo
        tmpInfo = excitInfo
        tmpInfo%fullEnd = excitInfo%firstEnd
        tmpInfo%currentGen = excitInfo%firstGen

        call calcAllExcitations(ilut, csf_i, tmpInfo, posSwitches, negSwitches, &
                                .false., tempExcits, tmpNum)

        ! then change tmpInfo to deal with second excitaion properly
        tmpInfo%fullStart = excitInfo%secondStart
        tmpInfo%fullEnd = excitInfo%fullEnd
        tmpInfo%currentGen = excitInfo%lastGen

        nExcits = 0
        ! have to allocate excitation list to worst casce
        nOpen1 = count_open_orbs_ij(csf_i, excitInfo%fullStart, excitInfo%firstEnd)
        nOpen2 = count_open_orbs_ij(csf_i, excitInfo%secondStart, excitInfo%fullEnd)

        nMax = 4 + 2**(nOpen1 + nOpen2 + 2)
        allocate(tmp_excitations(0:nifguga, nMax), stat=ierr)

        ! and for all created excitations i have to calc. all possible tops
        do iEx = 1, tmpNum
            call calcAllExcitations(tempExcits(:, iEx), csf_i, tmpInfo, posSwitches, &
                                    negSwitches, .false., tempExcits2, tmpNum2)

            ! have to reencode matrix element of tempexcits(:,iEx) as it is
            ! overwritten in the call allexcitation routine
            do iEx2 = 1, tmpNum2
                call update_matrix_element(tempExcits2(:, iEx2), &
                                           extract_matrix_element(tempExcits(:, iEx), 1), 1)

                ! also use the deltaB value of the finished excitations to
                ! handle IC in the NECI code correctly...
                call setDeltaB(2, tempExcits2(:, iEx2))
            end do

            ! and add it up to one list (maybe i can set sorted to true?
            ! since no repetitions in both...
            call add_guga_lists(nExcits, tmpNum2, tmp_excitations, tempExcits2)
        end do

        deallocate(tempExcits)
        deallocate(tempExcits2)

        j = 1
        do i = 1, nExcits
            if (near_zero(extract_matrix_element(tmp_excitations(:, i), 1))) cycle

            tmp_excitations(:, j) = tmp_excitations(:, i)

            j = j + 1
        end do

        nExcits = j - 1

        allocate(excitations(0:nifguga, nExcits), stat=ierr)
        excitations = tmp_excitations(:, 1:nExcits)

        deallocate(tmp_excitations)

    end subroutine calcNonOverlapDouble