findAlphaBetaOrbs Function

    function findAlphaBetaOrbs(symmax, nbasis)  result(idxAlphaBetaOrbs)
        use SymExcitDataMod, only: OrbClassCount
        use GenRandSymExcitNUMod, only: ClassCountInd
        integer, intent(in) :: symmax, nbasis
        integer :: idxAlphaBetaOrbs(nbasis)

        integer :: i, z, iSym, totEl, totOrbs
        integer :: alphaOrbs(symmax), betaOrbs(symmax), orbsSym(symmax+1)
        logical :: checkOccOrb

        orbsSym(:) = 0
        i = 1
        do iSym = 0, symmax-1
            i = i + 1
            orbsSym(i) = OrbClassCount(ClassCountInd(1, iSym, 0))*2 + orbsSym(i - 1)
        end do
        ! loop to find all the occupied alpha spin orbitals
        totEl = 0
        do iSym = 1, symmax
            alphaOrbs(iSym) = 0
            do z = 1, nel
                if (projEDet(z, 1) > orbsSym(iSym) .and. projEDet(z, 1) <= orbsSym(iSym + 1)) then
                    if (MOD(projEDet(z, 1), 2) == 1) then
                        alphaOrbs(iSym) = alphaOrbs(iSym) + 1
                        idxAlphaBetaOrbs(totEl + alphaOrbs(iSym)) = projEDet(z, 1)
                    end if
                end if
            end do
            totEl = totEl + alphaOrbs(iSym)
            ! loop to find all the occupied beta spin orbitals
            betaOrbs(iSym) = 0
            do z = 1, nel
                if (projEDet(z, 1) > orbsSym(iSym) .and. projEDet(z, 1) <= orbsSym(iSym + 1)) then
                    if (MOD(projEDet(z, 1), 2) == 0) then
                        betaOrbs(iSym) = betaOrbs(iSym) + 1
                        idxAlphaBetaOrbs(totEl + betaOrbs(iSym)) = projEDet(z, 1)
                    end if
                end if
            end do
            totEl = totEl + betaOrbs(iSym)
        end do
        if (totEl /= nel) write (stdout, *) 'WARNING: not matching number of electrons!'
        ! loop to find all the non-occupied alpha spin orbitals
        totOrbs = 0
        do iSym = 1, symmax
            alphaOrbs(iSym) = 0
            do i = orbsSym(iSym) + 1, orbsSym(iSym + 1)
                checkOccOrb = .false.
                do z = 1, nel
                    if (i == projEDet(z, 1)) checkOccOrb = .true.
                end do
                if (MOD(i, 2) == 1 .and. .not.checkOccOrb) then
                    alphaOrbs(iSym) = alphaOrbs(iSym) + 1
                    idxAlphaBetaOrbs(nel + totOrbs + alphaOrbs(iSym)) = i
                end if
            end do
            totOrbs = totOrbs + alphaOrbs(iSym)
            ! loop to find all the non-occupied beta spin orbitals
            betaOrbs(iSym) = 0
            do i = orbsSym(iSym) + 1, orbsSym(iSym + 1)
                checkOccOrb = .false.
                do z = 1, nel
                    if (i == projEDet(z, 1)) checkOccOrb = .true.
                end do
                if (MOD(i, 2) == 0 .and. .not.checkOccOrb) then
                    betaOrbs(iSym) = betaOrbs(iSym) + 1
                    idxAlphaBetaOrbs(nel + totOrbs + betaOrbs(iSym)) = i
                end if
            end do
            totOrbs = totOrbs + betaOrbs(iSym)
        end do
        if (nel + totOrbs /= nbasis) write (stdout, *) 'WARNING: not matching number of orbitals!'
    end function findAlphaBetaOrbs