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

evolution for a single Trotter step.

    Notes: follows Equation 9 of Poulin et al.'s work in "The Trotter Step
        Size Required for Accurate Quantum Simulation of Quantum Chemistry",
        applied to the "stagger"-based Trotter step for detailed in
        Kivlichan et al., "Quantum Simulation of Electronic Structure with
        Linear Depth and Connectivity", arxiv:1711.04789.
    """
    more_info = bool(indices)
    n_terms = len(terms)

    if verbose:
        import time
        start = time.time()

    error_operator = FermionOperator.zero()
    for beta in range(n_terms):
        if verbose and beta % (n_terms // 30) == 0:
            print('%4.3f percent done in' % (
                (float(beta) / n_terms) ** 3 * 100), time.time() - start)

        for alpha in range(beta + 1):
            for alpha_prime in range(beta):
                # If we have pre-computed info on indices, use it to determine
                # trivial double commutation.
                if more_info:
                    if (not
                        trivially_double_commutes_dual_basis_using_term_info(
                            indices[alpha], indices[beta],
                            indices[alpha_prime], is_hopping_operator[alpha],
                            is_hopping_operator[beta],
                            is_hopping_operator[alpha_prime], jellium_only)):

((left, 1), (left, 0)), 0.0))

            # Calculate the right number operator, right^ right.
            right_number_operator = FermionOperator(
                ((right, 1), (right, 0)), hamiltonian.terms.get(
                    ((right, 1), (right, 0)), 0.0))

            # Divide single-number terms by n_qubits-1 to avoid over-accounting
            # for the interspersed rotations. Each qubit is swapped n_qubits-1
            # times total.
            left_number_operator /= (n_qubits - 1)
            right_number_operator /= (n_qubits - 1)

        else:
            left_number_operator = FermionOperator.zero()
            right_number_operator = FermionOperator.zero()

        # If the overall hopping operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a hopping operator.
        if not (left_hopping_operator +
                right_hopping_operator) == FermionOperator.zero():
            terms_in_layer.append(left_hopping_operator +
                                  right_hopping_operator)
            indices_in_layer.append(set((left, right)))
            is_hopping_operator_in_layer.append(True)

        # If the overall number operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a number operator.
        if not (two_number_operator + left_number_operator +
                right_number_operator) == FermionOperator.zero():
            terms_in_layer.append(two_number_operator +
                                  left_number_operator +

def _fourier_transform_helper(hamiltonian,
                              grid,
                              spinless,
                              phase_factor,
                              vec_func_1,
                              vec_func_2):
    hamiltonian_t = FermionOperator.zero()
    normalize_factor = numpy.sqrt(1.0 / float(grid.num_points))

    for term in hamiltonian.terms:
        transformed_term = FermionOperator.identity()
        for ladder_op_mode, ladder_op_type in term:
            indices_1 = grid.grid_indices(ladder_op_mode, spinless)
            vec1 = vec_func_1(indices_1)
            new_basis = FermionOperator.zero()
            for indices_2 in grid.all_points_indices():
                vec2 = vec_func_2(indices_2)
                spin = None if spinless else ladder_op_mode % 2
                orbital = grid.orbital_id(indices_2, spin)
                exp_index = phase_factor * 1.0j * numpy.dot(vec1, vec2)
                if ladder_op_type == 1:
                    exp_index *= -1.0

((left, 1), (left, 0)), hamiltonian.terms.get(
                    ((left, 1), (left, 0)), 0.0))

            # Calculate the right number operator, right^ right.
            right_number_operator = FermionOperator(
                ((right, 1), (right, 0)), hamiltonian.terms.get(
                    ((right, 1), (right, 0)), 0.0))

            # Divide single-number terms by n_qubits-1 to avoid over-accounting
            # for the interspersed rotations. Each qubit is swapped n_qubits-1
            # times total.
            left_number_operator /= (n_qubits - 1)
            right_number_operator /= (n_qubits - 1)

        else:
            left_number_operator = FermionOperator.zero()
            right_number_operator = FermionOperator.zero()

        # If the overall hopping operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a hopping operator.
        if not (left_hopping_operator +
                right_hopping_operator) == FermionOperator.zero():
            terms_in_layer.append(left_hopping_operator +
                                  right_hopping_operator)
            indices_in_layer.append(set((left, right)))
            is_hopping_operator_in_layer.append(True)

        # If the overall number operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a number operator.
        if not (two_number_operator + left_number_operator +
                right_number_operator) == FermionOperator.zero():
            terms_in_layer.append(two_number_operator +

((right, 1), (right, 0)), 0.0))

            # Divide single-number terms by n_qubits-1 to avoid over-accounting
            # for the interspersed rotations. Each qubit is swapped n_qubits-1
            # times total.
            left_number_operator /= (n_qubits - 1)
            right_number_operator /= (n_qubits - 1)

        else:
            left_number_operator = FermionOperator.zero()
            right_number_operator = FermionOperator.zero()

        # If the overall hopping operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a hopping operator.
        if not (left_hopping_operator +
                right_hopping_operator) == FermionOperator.zero():
            terms_in_layer.append(left_hopping_operator +
                                  right_hopping_operator)
            indices_in_layer.append(set((left, right)))
            is_hopping_operator_in_layer.append(True)

        # If the overall number operator isn't close to zero, append it.
        # Include the indices it acts on and that it's a number operator.
        if not (two_number_operator + left_number_operator +
                right_number_operator) == FermionOperator.zero():
            terms_in_layer.append(two_number_operator +
                                  left_number_operator +
                                  right_number_operator)
            terms_in_layer[-1].compress()

            indices_in_layer.append(set((left, right)))
            is_hopping_operator_in_layer.append(False)

the commutator. All three operators must be of the same type.
        indices2, indices3 (set): The indices op2 and op3 act on.
        is_hopping_operator2 (bool): Whether op2 is a hopping operator.
        is_hopping_operator3 (bool): Whether op3 is a hopping operator.

    Returns:
        The double commutator of the given operators.
    """
    if is_hopping_operator2 and is_hopping_operator3:
        indices2 = set(indices2)
        indices3 = set(indices3)
        # Determine which indices both op2 and op3 act on.
        try:
            intersection, = indices2.intersection(indices3)
        except ValueError:
            return FermionOperator.zero()

        # Remove the intersection from the set of indices, since it will get
        # cancelled out in the final result.
        indices2.remove(intersection)
        indices3.remove(intersection)

        # Find the indices of the final output hopping operator.
        index2, = indices2
        index3, = indices3
        coeff2 = op2.terms[list(op2.terms)[0]]
        coeff3 = op3.terms[list(op3.terms)[0]]
        commutator23 = (
            FermionOperator(((index2, 1), (index3, 0)), coeff2 * coeff3) +
            FermionOperator(((index3, 1), (index2, 0)), -coeff2 * coeff3))
    else:
        commutator23 = normal_ordered(commutator(op2, op3))

def _fourier_transform_helper(hamiltonian,
                              grid,
                              spinless,
                              phase_factor,
                              vec_func_1,
                              vec_func_2):
    hamiltonian_t = FermionOperator.zero()
    normalize_factor = numpy.sqrt(1.0 / float(grid.num_points))

    for term in hamiltonian.terms:
        transformed_term = FermionOperator.identity()
        for ladder_op_mode, ladder_op_type in term:
            indices_1 = grid.grid_indices(ladder_op_mode, spinless)
            vec1 = vec_func_1(indices_1)
            new_basis = FermionOperator.zero()
            for indices_2 in grid.all_points_indices():
                vec2 = vec_func_2(indices_2)
                spin = None if spinless else ladder_op_mode % 2
                orbital = grid.orbital_id(indices_2, spin)
                exp_index = phase_factor * 1.0j * numpy.dot(vec1, vec2)
                if ladder_op_type == 1:
                    exp_index *= -1.0

                element = FermionOperator(((orbital, ladder_op_type),),
                                          numpy.exp(exp_index))
                new_basis += element

            new_basis *= normalize_factor
            transformed_term *= new_basis

        # Coefficient.