CalcPGenHPHF Subroutine

    subroutine CalcPGenHPHF(nI, iLutnI, nJ, iLutnJ, ex, ClassCount, ClassCountUnocc, pDoubles, pGen, tSameFunc)
        integer, intent(in) :: nI(nel)
        integer(kind=n_int), intent(in) :: iLutnI(0:niftot), iLutnJ(0:niftot)
        integer, intent(in) :: ClassCount(ScratchSize), ClassCountUnocc(ScratchSize)
        integer, intent(in) :: nJ(nel), ex(2, maxExcit)
        real(dp), intent(in) :: pDoubles
        real(dp), intent(out) :: pGen
        logical, intent(out) :: tSameFunc
        logical :: tSign, tSwapped
        real(dp) :: pGen2
        integer :: ic
        integer :: Ex2(2, maxExcit), nJ_loc(nel), nJ2(nel)
        integer(kind=n_int) :: iLutnJ_loc(0:niftot), iLutnJ2(0:niftot)
#ifdef DEBUG_
        character(*), parameter :: this_routine = "CalcPGenHPHF"
#endif
        tSameFunc = .false.
        pGen = 0.0_dp

        IF(TestClosedShellDet(iLutnJ)) THEN
            !nJ is CS, therefore, only one way of generating it.
            ic = FindBitExcitLevel(iLutnI, iLutnJ, 2)
            if(ic == 0) then
                tSameFunc = .true.
                return
            end if
            if(ic <= max_ex_level) then
                call CalcNonUniPGen(nI, ilutnI, ex, ic, ClassCount, ClassCountUnocc, pDoubles, pGen)
            end if
        else
            !nJ is openshell. Add the probabilities of generating each pair (if both connected)
            nJ_loc = nJ
            iLutnJ_loc = iLutnJ
            CALL ReturnAlphaOpenDet(nJ_loc, nJ2, iLutnJ_loc, iLutnJ2, .true., .true., tSwapped)

            !First find nI -> nJ
            ic = FindBitExcitLevel(iLutnI, iLutnJ_loc, 2)
            if(ic == 0) then
                tSameFunc = .true.
                return
            end if
            if(ic <= max_ex_level) then
                if(.not. tSwapped) then
                    !ex is correct for this excitation
                    call CalcNonUnipGen(nI, ilutnI, ex, ic, ClassCount, ClassCountUnocc, pDoubles, pGen)
                else
                    Ex2(1, 1) = ic
                    call GetBitExcitation(iLutnI, iLutnJ_loc, Ex2, tSign)
                    call CalcNonUnipGen(nI, ilutnI, Ex2, ic, ClassCount, ClassCountUnocc, pDoubles, pGen)
                end if
            end if

            !Now consider nI -> nJ2 and add the probabilities
            ic = FindBitExcitLevel(iLutnI, iLutnJ2, 2)
            if(ic == 0) then
                tSameFunc = .true.
                return
            end if
            if(ic <= max_ex_level) then
                if(tSwapped) then
                    !ex is correct for this excitation
                    call CalcNonUnipGen(nI, ilutnI, ex, ic, ClassCount, ClassCountUnocc, pDoubles, pGen2)
                else
                    Ex2(1, 1) = ic
                    call GetBitExcitation(iLutnI, iLutnJ2, Ex2, tSign)
                    call CalcNonUnipGen(nI, ilutnI, Ex2, ic, ClassCount, ClassCountUnocc, pDoubles, pGen2)
                end if
                pGen = pGen + pGen2
            end if
        end if

    end subroutine CalcPGenHPHF