import mpmath


def energy_el(P, F, H):
    """
    Compute electronic energy.

    INPUT:
        P: Density matrix
        F: Fock matrix
        H: Core Hamiltonian

    Source:
        Modern Quantum Chemistry
        Szabo and Ostlund
        Dover
        1989
    """

    # Size of the basis set
    K = P.shape[0]

    E = 0

    for i in range(K):
        for j in range(K):
            E += mpmath.mp.mpf("0.5") * P[i, j] * (H[i, j] + F[i, j])

    return E