hphf_diag_helement Function

    function hphf_diag_helement(nI, iLutnI) result(hel)

        ! Find the diagonal HElment for a half-projected hartree-fock
        ! determinant.
        !
        ! In:  nI      - Determinant to consider
        !      iLutnI  - Bit representation of I
        ! Ret: hel   - The calculated matrix element

        integer, intent(in) :: nI(nel)
        integer(kind=n_int), intent(in) :: iLutnI(0:NIfTot)
        HElement_t(dp) :: hel

        integer(kind=n_int) :: iLutnI2(0:NIfTot)
        integer :: ExcitLevel, OpenOrbs
        HElement_t(dp) :: MatEl2
        integer :: nJ(nel)

        if(t_lattice_model) then
            hel = get_helement_lattice(nI, nI)
        else
            hel = sltcnd_excit(nI, Excite_0_t())
        end if

        if(.not. TestClosedShellDet(iLutnI)) then
            ! <i|H|i> = <j|H|j>, so no need to calculate both.
            ! <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. In the case of the
            ! antisymmetric pair, the cross term is subtracted.

            ! See if there is a cross-term
            call FindExcitBitDetSym(iLutnI, iLutnI2)
            ExcitLevel = FindBitExcitLevel(iLutnI, iLutnI2, 2)
            if(ExcitLevel <= max_ex_level) then
                call CalcOpenOrbs(iLutnI, OpenOrbs)
                if(t_lattice_model) then
                    call decode_bit_det(nJ, iLutnI2)
                    ! here i am really not sure about hermiticity..
                    MatEl2 = get_helement_lattice(nI, nJ)
                    ! do i need a hermitian version of that here?
                else
                    MatEl2 = sltcnd(nI, iLutnI, iLutnI2)
                end if

                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

        hel = hel + (ECore)
    end function hphf_diag_helement