How to use the openfermion.FermionOperator.zero function in openfermion

for. Identity is arbitrarily chosen
                                       to be part of the potential.

    Returns:
        Tuple of (potential_terms, kinetic_terms). Both elements of the tuple
        are numpy arrays of FermionOperators.
    """
    if not isinstance(hamiltonian, FermionOperator):
        try:
            hamiltonian = normal_ordered(get_fermion_operator(hamiltonian))
        except TypeError:
            raise TypeError('hamiltonian must be either a FermionOperator '
                            'or DiagonalCoulombHamiltonian.')

    potential = FermionOperator.zero()
    kinetic = FermionOperator.zero()

    for term, coeff in iteritems(hamiltonian.terms):
        acted = set(term[i][0] for i in range(len(term)))
        if len(acted) == len(term) / 2:
            potential += FermionOperator(term, coeff)
        else:
            kinetic += FermionOperator(term, coeff)

    potential_terms = numpy.array(
        [FermionOperator(term, coeff)
         for term, coeff in iteritems(potential.terms)])

    kinetic_terms = numpy.array(
        [FermionOperator(term, coeff)
         for term, coeff in iteritems(kinetic.terms)])

Args:
        side_length: The side length of the 2D jellium grid. There are
            side_length ** 2 qubits, and O(side_length ** 4) terms in the
            Hamiltonian.

    Returns:
        runtime_commutator: The time it takes to compute a commutator, after
            partitioning the terms and normal ordering, using the regular
            commutator function.
        runtime_diagonal_commutator: The time it takes to compute the same
            commutator using methods restricted to diagonal Coulomb operators.
    """
    hamiltonian = normal_ordered(jellium_model(Grid(2, side_length, 1.),
                                               plane_wave=False))

    part_a = FermionOperator.zero()
    part_b = FermionOperator.zero()
    add_to_a_or_b = 0  # add to a if 0; add to b if 1
    for term, coeff in iteritems(hamiltonian.terms):
        # Partition terms in the Hamiltonian into part_a or part_b
        if add_to_a_or_b:
            part_a += FermionOperator(term, coeff)
        else:
            part_b += FermionOperator(term, coeff)
        add_to_a_or_b ^= 1

    start = time.time()
    _ = normal_ordered(commutator(part_a, part_b))
    end = time.time()
    runtime_commutator = end - start

    start = time.time()

Notes:
        Follows Eq 9 of Poulin et al., arXiv:1406.4920, applied to the
        Trotter step detailed in Kivlichan et al., arxiv:1711.04789.
    """
    single_terms = numpy.array(
        simulation_ordered_grouped_low_depth_terms_with_info(
            hamiltonian,
            external_potential_at_end=external_potential_at_end)[0])

    # Cache the halved terms for use in the second commutator.
    halved_single_terms = single_terms / 2.0

    term_mode_mask = bit_mask_of_modes_acted_on_by_fermionic_terms(
        single_terms, count_qubits(hamiltonian))

    error_operator = FermionOperator.zero()

    for beta, term_beta in enumerate(single_terms):
        modes_acted_on_by_term_beta = set()
        for beta_action in term_beta.terms:
            modes_acted_on_by_term_beta.update(
                set(operator[0] for operator in beta_action))

        beta_mode_mask = numpy.logical_or.reduce(
            [term_mode_mask[mode] for mode in modes_acted_on_by_term_beta])

        # alpha_prime indices that could have a nonzero commutator, i.e.
        # there's overlap between the modes the corresponding terms act on.
        valid_alpha_primes = numpy.where(beta_mode_mask)[0]

        # Only alpha_prime < beta enters the error operator; filter for this.
        valid_alpha_primes = valid_alpha_primes[valid_alpha_primes < beta]

side_length: The side length of the 2D jellium grid. There are
            side_length ** 2 qubits, and O(side_length ** 4) terms in the
            Hamiltonian.

    Returns:
        runtime_commutator: The time it takes to compute a commutator, after
            partitioning the terms and normal ordering, using the regular
            commutator function.
        runtime_diagonal_commutator: The time it takes to compute the same
            commutator using methods restricted to diagonal Coulomb operators.
    """
    hamiltonian = normal_ordered(jellium_model(Grid(2, side_length, 1.),
                                               plane_wave=False))

    part_a = FermionOperator.zero()
    part_b = FermionOperator.zero()
    add_to_a_or_b = 0  # add to a if 0; add to b if 1
    for term, coeff in iteritems(hamiltonian.terms):
        # Partition terms in the Hamiltonian into part_a or part_b
        if add_to_a_or_b:
            part_a += FermionOperator(term, coeff)
        else:
            part_b += FermionOperator(term, coeff)
        add_to_a_or_b ^= 1

    start = time.time()
    _ = normal_ordered(commutator(part_a, part_b))
    end = time.time()
    runtime_commutator = end - start

    start = time.time()
    _ = commutator_ordered_diagonal_coulomb_with_two_body_operator(

separate the potential and kinetic terms
                                       for. Identity is arbitrarily chosen
                                       to be part of the potential.

    Returns:
        Tuple of (potential_terms, kinetic_terms). Both elements of the tuple
        are numpy arrays of FermionOperators.
    """
    if not isinstance(hamiltonian, FermionOperator):
        try:
            hamiltonian = normal_ordered(get_fermion_operator(hamiltonian))
        except TypeError:
            raise TypeError('hamiltonian must be either a FermionOperator '
                            'or DiagonalCoulombHamiltonian.')

    potential = FermionOperator.zero()
    kinetic = FermionOperator.zero()

    for term, coeff in iteritems(hamiltonian.terms):
        acted = set(term[i][0] for i in range(len(term)))
        if len(acted) == len(term) / 2:
            potential += FermionOperator(term, coeff)
        else:
            kinetic += FermionOperator(term, coeff)

    potential_terms = numpy.array(
        [FermionOperator(term, coeff)
         for term, coeff in iteritems(potential.terms)])

    kinetic_terms = numpy.array(
        [FermionOperator(term, coeff)
         for term, coeff in iteritems(kinetic.terms)])

# Cache halved potential and kinetic terms for the second commutator.
    halved_potential_terms = potential_terms / 2.0
    halved_kinetic_terms = kinetic_terms / 2.0

    # Assign the outer term of the second commutator based on the ordering.
    outer_potential_terms = (halved_potential_terms if order == 'T+V' else
                             potential_terms)
    outer_kinetic_terms = (halved_kinetic_terms if order == 'V+T' else
                           kinetic_terms)

    potential_mask = bit_mask_of_modes_acted_on_by_fermionic_terms(
        potential_terms, n_qubits)
    kinetic_mask = bit_mask_of_modes_acted_on_by_fermionic_terms(
        kinetic_terms, n_qubits)

    error_operator = FermionOperator.zero()

    for potential_term in potential_terms:
        modes_acted_on_by_potential_term = set()

        for potential_term_action in potential_term.terms:
            modes_acted_on_by_potential_term.update(
                set(operator[0] for operator in potential_term_action))

        if not modes_acted_on_by_potential_term:
            continue

        potential_term_mode_mask = numpy.logical_or.reduce(
            [kinetic_mask[mode] for mode in modes_acted_on_by_potential_term])

        for kinetic_term in kinetic_terms[potential_term_mode_mask]:
            inner_commutator_term = (