get_helement_normal Function

    function get_helement_normal(nI, nJ, iLutI, iLutJ, ICret) result(hel)

        ! Get the matrix element of the hamiltonian.
        !
        ! In:  nI, nJ       - The determinants to consider
        !      iLutI, iLutJ - Bit representations of I,J (optional, helpful)
        ! Out: ICret        - The number of orbitals I,J differ by
        ! Ret: hel          - The desired matrix element.

        integer, intent(in) :: nI(nel), nJ(nel)
        integer(kind=n_int), intent(in), optional :: iLutI(0:NIfTot), iLutJ(0:NIfTot)
        integer, intent(out), optional :: ICret
        HElement_t(dp) :: hel

        character(*), parameter :: this_routine = 'get_helement_normal'
        integer :: ex(2, 2), IC
        integer(kind=n_int) :: ilut(0:NIfTot, 2)

        if (tGUGA) then

            if (all(nI == nJ)) then
                hel = calcDiagMatEleGUGA_nI(nI)
            else
                call stop_all(this_routine, "TODO: refactor guga matrix elements!")
            end if
            return
        end if

        if (tHPHFInts) &
            call stop_all(this_routine, "Should not be calling HPHF &
                          &integrals from here.")

        if (t_lattice_model) then
            if (present(ICret)) then
                ic = -1
                hel = get_helement_lattice(nI, nJ, ic)
                ICret = ic
            else
                hel = get_helement_lattice(nI, nJ)
            end if
            return
        end if

        if (tStoreAsExcitations .and. nI(1) == -1 .and. nJ(1) == -1) then

            ex(1, :) = nJ(4:5)
            ex(2, :) = nJ(6:7)
            hel = sltcnd_excit(nI, Excite_2_t(ex), .false.)
        end if

        if (present(iLutJ)) then
            hel = sltcnd(nI, iLutI, iLutJ, IC)
        else
            call EncodeBitDet(nI, iLut(:, 1))
            call EncodeBitdet(nJ, iLut(:, 2))
            ! TODO: This is not an ideal place to end up...
            hel = sltcnd(nI, iLut(:, 1), ilut(:, 2), IC)
        end if

        ! Add in ECore for a diagonal element
        if (IC == 0) then
            hel = hel + (ECore)
        else if (tStoquastize) then
            hel = -abs(hel)
        end if

        ! If requested, return IC
        if (present(ICret)) then
            ICret = IC
        end if

    end function get_helement_normal