SumFock Function

    pure function SumFock(nI, HFDet) result(hel)

        ! This just calculates the sum of the Fock energies
        ! by considering the one-electron integrals and
        ! the double-counting contribution
        ! to the diagonal matrix elements. This is subtracted from
        ! the sum of the fock energies to calculate diagonal
        ! matrix elements, or added to the sum of the 1-electron
        ! integrals. The HF determinant needs to be supplied.

        integer, intent(in) :: nI(nel), HFDet(nel)
        HElement_t(dp) :: hel
        integer :: i

        hel = h_cast(0._dp)
        do i = 1, nEl
            hel = hel + CalcFockOrbEnergy(nI(i), HFDet)
        end do
    end function SumFock