isDouble Function

    pure function isDouble(nI, iOrb) result(flag)
        ! returns a logical .true. if spinorbital iOrb is doubly occupied
        ! and .false. elsewise.
        integer, intent(in) :: nI(nEl), iOrb
        logical :: flag

        integer :: pair

        flag = .false.
        pair = ab_pair(nI(iOrb))

        ! ask simon about the ordering in nI. If its not always ordered I
        ! have to to a quicksearch for pair in nI. If its ordered I can
        ! just check adjacent entries in nI

        ! assume ordered for now!
        if (iOrb == 1) then
            ! special conditions if iOrb is 1
            flag = (nI(2) == pair)
        else if (iOrb == nEl) then
            flag = (nI(nEl - 1) == pair)
        else
            flag = ((nI(iOrb - 1) == pair) .or. (nI(iOrb + 1) == pair))
        end if

    end function isDouble