correct_hdiag_hphf Subroutine

    subroutine correct_hdiag_hphf(nJ, iLutnJ, 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) :: nJ(nel)
        integer(n_int), intent(in) :: iLutnJ(0:NIfTot)
        HElement_t(dp), intent(inout) :: hel

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

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

                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 subroutine correct_hdiag_hphf