CalcFockOrbEnergy Function

    pure function CalcFockOrbEnergy(Orb, HFDet) result(hel)
        ! This calculates the orbital fock energy from
        ! the one- and two- electron integrals. This
        ! requires a knowledge of the HF determinant.
        !In: Orbital (Spin orbital notation)
        !In: HFDet (HF Determinant)
        integer, intent(in) :: HFDet(nel), Orb
        integer :: idHF(NEl), idOrb, j
        HElement_t(dp) :: hel

        ! Obtain the spatial rather than spin indices if required
        idOrb = gtID(Orb)
        idHF = gtID(HFDet)

        !GetTMATEl works with spin orbitals
        hel = GetTMATEl(Orb, Orb)
        ! Sum in the two electron contributions.
        do j = 1, nel
            hel = hel + get_umat_el(idOrb, idHF(j), idOrb, idHF(j))
        end do

        ! Exchange contribution only considered if tExch set.
        ! This is only separated from the above loop to keep "if (tExch)" out
        ! of the tight loop for efficiency.
        do j = 1, nel
            ! Exchange contribution is zero if I,J are alpha/beta
            if (tReltvy .or. (G1(Orb)%Ms == G1(HFDet(j))%Ms)) then
                hel = hel - get_umat_el(idOrb, idHF(j), idHF(j), idOrb)
            end if
        end do

    end function CalcFockOrbEnergy