get_hdiag_bare_hphf Function

    function get_hdiag_bare_hphf(nI, iLutnI, hel_old) result(hel)

        ! <X|H|X> = 1/2 [ <i|H|i> + <j|H|j> ] + <i|H|j> where i and j are
        ! the two spin-coupled dets which make up X.

        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: iLutnI(0:NIfTot)
        real(dp), intent(in) :: hel_old

        HElement_t(dp) :: hel

        integer(n_int) :: iLutnI2(0:NIfTot)
        integer :: ExcitLevel, OpenOrbs
        HElement_t(dp) :: MatEl2

        hel = hel_old

        if (.not. TestClosedShellDet(iLutnI)) then
            ! See if there is a cross-term. If there is, then remove this
            ! from hel_old to get the desired value
            call FindExcitBitDetSym(iLutnI, iLutnI2)
            ExcitLevel = FindBitExcitLevel(iLutnI, iLutnI2, 2)
            if (ExcitLevel <= 2) then
                call CalcOpenOrbs(iLutnI, OpenOrbs)
                MatEl2 = sltcnd(nI, iLutnI, iLutnI2)

                if (tOddS_HPHF) then
                    if (mod(OpenOrbs, 2) == 1) then
                        ! Subtract cross terms if determinant is antisymmetric.
                        hel = hel - MatEl2
                    else
                        hel = hel + MatEl2
                    end if
                else
                    if (mod(OpenOrbs, 2) == 1) then
                        ! Subtract cross terms if determinant is antisymmetric.
                        hel = hel + MatEl2
                    else
                        hel = hel - MatEl2
                    end if
                end if
            end if
        end if

    end function get_hdiag_bare_hphf