How to use the openfermion.utils.normal_ordered function in openfermion

Args:
        hamiltonian (FermionOperator): The Hamiltonian.

    Returns:
        A list of terms from the Hamiltonian in simulated order.

    Notes:
        Assumes the Hamiltonian is in the form discussed in Kivlichan
        et al., "Quantum Simulation of Electronic Structure with Linear
        Depth and Connectivity", arxiv:1711.04789. This constrains the
        Hamiltonian to have only hopping terms (i^ j + j^ i) and potential
        terms which are products of at most two number operators (n_i or
        n_i n_j).
    """
    n_qubits = count_qubits(hamiltonian)
    hamiltonian = normal_ordered(hamiltonian)
    terms = []

    for operators, coefficient in iteritems(hamiltonian.terms):
        terms += [FermionOperator(operators, coefficient)]

    return terms

first_round_swaps, second_round_swaps = ffft_swap_networks(system_size,
                                                               n_dimensions=2)

    # Create the full list of swaps.
    swap_mode_list = first_round_swaps + second_round_swaps
    swaps = [swap[0] for swap_layer in swap_mode_list for swap in swap_layer]

    # Swap adjacent modes for all swaps in both rounds.
    for mode in swaps:
        fermion_operator = swap_adjacent_fermionic_modes(
            fermion_operator, mode)
        fermion_operator = normal_ordered(fermion_operator)

    # Apply the Fourier transform to the reordered operator.
    fft_result = fourier_transform(fermion_operator, grid, spinless=True)
    fft_result = normal_ordered(fft_result)

    # Undo the second round of swaps to restore the FFT's bit-reversed order.
    swap_mode_list = second_round_swaps
    swaps = [swap[0] for swap_layer in swap_mode_list for swap in swap_layer]
    for mode in swaps:
        fft_result = swap_adjacent_fermionic_modes(fft_result, mode)
        fft_result = normal_ordered(fft_result)

    return fft_result

hartree_fock_state_index = numpy.sum(2 ** occupied_states)
        hartree_fock_state = numpy.zeros(2 ** count_qubits(hamiltonian),
                                         dtype=complex)
        hartree_fock_state[hartree_fock_state_index] = 1.0

    else:
        # Inverse Fourier transform the creation operators for the state to get
        # to the dual basis state, then use that to get the dual basis state.
        hartree_fock_state_creation_operator = FermionOperator.identity()
        for state in occupied_states[::-1]:
            hartree_fock_state_creation_operator *= (
                FermionOperator(((int(state), 1),)))
        dual_basis_hf_creation_operator = inverse_fourier_transform(
            hartree_fock_state_creation_operator, grid, spinless)

        dual_basis_hf_creation = normal_ordered(
            dual_basis_hf_creation_operator)

        # Initialize the HF state.
        hartree_fock_state = numpy.zeros(2 ** count_qubits(hamiltonian),
                                         dtype=complex)

        # Populate the elements of the HF state in the dual basis.
        for term in dual_basis_hf_creation.terms:
            index = 0
            for operator in term:
                index += 2 ** operator[0]
            hartree_fock_state[index] = dual_basis_hf_creation.terms[term]

    return hartree_fock_state

of this method when it encounters terms that are not represented
            by the DiagonalCoulombHamiltonian class, namely, terms that are
            not quadratic and not quartic of the form
            a^\dagger_p a_p a^\dagger_q a_q. If set to True, this method will
            simply ignore those terms. If False, then this method will raise
            an error if it encounters such a term. The default setting is False.
    """
    if not isinstance(fermion_operator, FermionOperator):
        raise TypeError('Input must be a FermionOperator.')

    if n_qubits is None:
        n_qubits = count_qubits(fermion_operator)
    if n_qubits < count_qubits(fermion_operator):
        raise ValueError('Invalid number of qubits specified.')

    fermion_operator = normal_ordered(fermion_operator)
    constant = 0.
    one_body = numpy.zeros((n_qubits, n_qubits), complex)
    two_body = numpy.zeros((n_qubits, n_qubits), float)

    for term, coefficient in fermion_operator.terms.items():
        # Ignore this term if the coefficient is zero
        if abs(coefficient) < EQ_TOLERANCE:
            continue

        if len(term) == 0:
            constant = coefficient
        else:
            actions = [operator[1] for operator in term]
            if actions == [1, 0]:
                p, q = [operator[0] for operator in term]
                one_body[p, q] = coefficient

hubbard_hamiltonian (FermionOperator): The Hamiltonian.
        original_ordering (list): The initial Jordan-Wigner canonical order.

    Returns:
        A 3-tuple of terms from the Hubbard Hamiltonian in order of
        simulation, the indices they act on, and whether they are hopping
        operators (both also in the same order).

    Notes:
        Assumes that the Hubbard model has spin and is on a 2D square
        aperiodic lattice. Uses the "stagger"-based Trotter step for the
        Hubbard model detailed in Kivlichan et al., "Quantum Simulation
        of Electronic Structure with Linear Depth and Connectivity",
        arxiv:1711.04789.
    """
    hamiltonian = normal_ordered(hubbard_hamiltonian)
    n_qubits = count_qubits(hamiltonian)
    side_length = int(numpy.sqrt(n_qubits / 2.0))

    ordered_terms = []
    ordered_indices = []
    ordered_is_hopping_operator = []

    original_ordering = list(range(n_qubits))
    for i in range(0, n_qubits - side_length, 2 * side_length):
        for j in range(2 * bool(i % (4 * side_length)), 2 * side_length, 4):
            original_ordering[i+j], original_ordering[i+j+1] = (
                original_ordering[i+j+1], original_ordering[i+j])

    input_ordering = list(original_ordering)

    # Follow odd-even transposition sort. In alternating steps, swap each even

Warning:
        Even assuming that each creation or annihilation operator appears
        at most a constant number of times in the original operator, the
        runtime of this method is exponential in the number of qubits.
    """
    if not isinstance(fermion_operator, FermionOperator):
        raise TypeError('Input must be a FermionOperator.')

    if n_qubits is None:
        n_qubits = count_qubits(fermion_operator)
    if n_qubits < count_qubits(fermion_operator):
        raise ValueError('Invalid number of qubits specified.')

    # Normal order the terms and initialize.
    fermion_operator = normal_ordered(fermion_operator)
    constant = 0.
    combined_hermitian_part = numpy.zeros((n_qubits, n_qubits), complex)
    antisymmetric_part = numpy.zeros((n_qubits, n_qubits), complex)

    # Loop through terms and assign to matrix.
    for term in fermion_operator.terms:
        coefficient = fermion_operator.terms[term]
        # Ignore this term if the coefficient is zero
        if abs(coefficient) < EQ_TOLERANCE:
            continue

        if len(term) == 0:
            # Constant term
            constant = coefficient
        elif len(term) == 2:
            ladder_type = [operator[1] for operator in term]

n_electrons (int): Number of electrons in the system.
        spinless (bool): Whether to use the spinless model or not.
        plane_wave (bool): Whether to return the Hartree-Fock state in
                           the plane wave (True) or dual basis (False).

    Notes:
        The jellium model is built up by filling the lowest-energy
        single-particle states in the plane-wave Hamiltonian until
        n_electrons states are filled.
    """
    from openfermion.hamiltonians import plane_wave_kinetic
    # Get the jellium Hamiltonian in the plane wave basis.
    # For determining the Hartree-Fock state in the PW basis, only the
    # kinetic energy terms matter.
    hamiltonian = plane_wave_kinetic(grid, spinless=spinless)
    hamiltonian = normal_ordered(hamiltonian)
    hamiltonian.compress()

    # The number of occupied single-particle states is the number of electrons.
    # Those states with the lowest single-particle energies are occupied first.
    occupied_states = lowest_single_particle_energy_states(
        hamiltonian, n_electrons)
    occupied_states = numpy.array(occupied_states)

    if plane_wave:
        # In the plane wave basis the HF state is a single determinant.
        hartree_fock_state_index = numpy.sum(2 ** occupied_states)
        hartree_fock_state = numpy.zeros(2 ** count_qubits(hamiltonian),
                                         dtype=complex)
        hartree_fock_state[hartree_fock_state_index] = 1.0

    else:

Warning:
        Even assuming that each creation or annihilation operator appears
        at most a constant number of times in the original operator, the
        runtime of this method is exponential in the number of qubits.
    """
    if not isinstance(fermion_operator, FermionOperator):
        raise TypeError('Input must be a FermionOperator.')

    if n_qubits is None:
        n_qubits = count_qubits(fermion_operator)
    if n_qubits < count_qubits(fermion_operator):
        raise ValueError('Invalid number of qubits specified.')

    # Normal order the terms and initialize.
    fermion_operator = normal_ordered(fermion_operator)
    constant = 0.
    one_body = numpy.zeros((n_qubits, n_qubits), complex)
    two_body = numpy.zeros((n_qubits, n_qubits,
                            n_qubits, n_qubits), complex)

    # Loop through terms and assign to matrix.
    for term in fermion_operator.terms:
        coefficient = fermion_operator.terms[term]
        # Ignore this term if the coefficient is zero
        if abs(coefficient) < EQ_TOLERANCE:
            continue

        # Handle constant shift.
        if len(term) == 0:
            constant = coefficient