add_core_en Subroutine

    subroutine add_core_en(matel, indices)
        HElement_t(dp), intent(inout) :: matel
        integer(int64), intent(in) :: indices(num_inds)

        integer :: counts
        integer(int64) :: idx(num_ex), ct
        real(dp) :: prefactor, delta
        integer(MPIArg) :: ierr

        if(t_freeze(indices)) then
            ! Count the number of different indices appearing
            counts = count_frozen_inds(indices)
            ! How many pairs of frozen orbitals do we have
            select case(counts)
            ! 1 => double excitation
            case(1)
                idx = frozen_double_entry(indices, prefactor)
                ! There are up to eight potential double excitations with four given indices
                ! to which the LMat entry contributes - this is because LMat is hermitian
                ! but UMat is not
                do ct = 1, num_ex
                    ! If the index is assigned, there is a duplicate frozen orb
                    if(idx(ct) > 0) then
                        ! Absorb the matrix element into UMat (with the corresponding prefactor
                        ! according to permutation and number of possible spin terms
                        ! UMat(idx(ct)) = UMat(idx(ct)) + prefactor*matel
                        delta = prefactor * matel
                        ! Use a remote update of UMat
                        call MPI_Accumulate(delta, 1, MPI_DOUBLE_PRECISION, 0, &
                                            idx(ct) - 1, 1, MPI_DOUBLE_PRECISION, MPI_SUM, umat_win, ierr)
                    endif
                end do
                ! 2 => single excitation
            case(2)
                idx(1:2) = frozen_single_entry(indices, prefactor)

                if(idx(1) > 0) then
                    ! Absorb the matrix element into TMat
                    ! TMat2D is indexed with spin orbs
                    call add_to_tmat(spatToSpinAlpha(idx(1)), spatToSpinAlpha(idx(2)))
                    call add_to_tmat(spatToSpinBeta(idx(1)), spatToSpinBeta(idx(2)))
                endif
                ! 3 => diagonal element
            case(3)
                call frozen_diagonal_entry(indices, prefactor)
                ECore_local = ECore_local + prefactor * matel
            end select
            ! Zero the matrix element for further usage (i.e. will not turn up anymore)
            matel = 0.0_dp
        endif

    contains

        subroutine add_to_tmat(ind1, ind2)
            integer(int64), intent(in) :: ind1, ind2

            TMat_local(ind1, ind2) = TMat_local(ind1, ind2) + prefactor * matel
            ! TMat is hermitian, as is the 3-body term
            if(ind1 /= ind2) &
                TMat_local(ind2, ind1) = TMat_local(ind2, ind1) + prefactor * matel
        end subroutine add_to_tmat

    end subroutine add_core_en