How to use the pyquil.gates.RZ function in pyquil

def test_larger_qaoa_circuit(qvm):
    square_qaoa_circuit = [Hgate(0), Hgate(1), Hgate(2), Hgate(3),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           CNOTgate(0, 1),
                           RZgate(0.78564882604980579)(1),
                           CNOTgate(0, 1),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           CNOTgate(0, 3),
                           RZgate(0.78564882604980579)(3),
                           CNOTgate(0, 3),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           Xgate(0),
                           PHASEgate(0.3928244130249029)(0),
                           CNOTgate(1, 2),
                           RZgate(0.78564882604980579)(2),
                           CNOTgate(1, 2),
                           Xgate(0),

List of the single qubits to apply rotations on.

    Returns
    -------
    Program
        Parametric Quil code containing the Z-Rotations.

    """
    p = Program()

    if len(qubit_pairs) != gammas_pairs.declared_size:
        raise ValueError("zz_rotation_angles must have the same length as qubits_pairs")

    for gamma_pair, qubit_pair in zip(gammas_pairs, qubit_pairs):
        p.inst(RZ(2 * gamma_pair, qubit_pair[0]))
        p.inst(RZ(2 * gamma_pair, qubit_pair[1]))
        p.inst(CPHASE(-4 * gamma_pair, qubit_pair[0], qubit_pair[1]))

    if gammas_singles.declared_size != len(qubit_singles):
        raise ValueError("z_rotation_angles must have the same length as qubit_singles")

    for gamma_single, qubit in zip(gammas_singles, qubit_singles):
        p.inst(RZ(2 * gamma_single, qubit))

    return p

"""
    p = Program()

    if len(qubit_pairs) != gammas_pairs.declared_size:
        raise ValueError("zz_rotation_angles must have the same length as qubits_pairs")

    for gamma_pair, qubit_pair in zip(gammas_pairs, qubit_pairs):
        p.inst(RZ(2 * gamma_pair, qubit_pair[0]))
        p.inst(RZ(2 * gamma_pair, qubit_pair[1]))
        p.inst(CPHASE(-4 * gamma_pair, qubit_pair[0], qubit_pair[1]))

    if gammas_singles.declared_size != len(qubit_singles):
        raise ValueError("z_rotation_angles must have the same length as qubit_singles")

    for gamma_single, qubit in zip(gammas_singles, qubit_singles):
        p.inst(RZ(2 * gamma_single, qubit))

    return p

change_to_z_basis.inst(RX(np.pi / 2.0, index))
            change_to_original_basis.inst(RX(-np.pi / 2.0, index))

        elif "I" == op:
            continue

        if prev_index is not None:
            cnot_seq.inst(CNOT(prev_index, index))

        prev_index = index
        highest_target_index = index

    # building rotation circuit
    quil_prog += change_to_z_basis
    quil_prog += cnot_seq
    quil_prog.inst(RZ(2.0 * pauli_term.coefficient * param, highest_target_index))
    quil_prog += reverse_hack(cnot_seq)
    quil_prog += change_to_original_basis

    return quil_prog

def _one_q_sic_prep(index, qubit):
    """Prepare the index-th SIC basis state."""
    if index == 0:
        return Program()

    theta = 2 * np.arccos(1 / np.sqrt(3))
    zx_plane_rotation = Program([RX(-pi / 2, qubit), RZ(theta - pi, qubit), RX(-pi / 2, qubit)])

    if index == 1:
        return zx_plane_rotation

    elif index == 2:
        return zx_plane_rotation + RZ(-2 * pi / 3, qubit)

    elif index == 3:
        return zx_plane_rotation + RZ(2 * pi / 3, qubit)

    raise ValueError(f"Bad SIC index: {index}")

:param qubits: A list of ints corresponding to the qubits to operate on. The operator
                   operates on bistrings of the form |qubits[0], ..., qubits[-1]>.
    """
    p = pq.Program()

    if len(qubits) == 1:
        p.inst(H(qubits[0]))
        p.inst(Z(qubits[0]))
        p.inst(H(qubits[0]))

    else:
        p.inst(map(X, qubits))
        p.inst(H(qubits[-1]))
        p.inst(RZ(-np.pi)(qubits[0]))
        p += n_qubit_control(qubits[:-1], qubits[-1], np.array([[0, 1], [1, 0]]), "NOT")
        p.inst(RZ(-np.pi)(qubits[0]))
        p.inst(H(qubits[-1]))
        p.inst(map(X, qubits))
    return p

qubit_singles:
        List of the single qubits to apply rotations on.

    Returns
    -------
    Program
        Parametric Quil code containing the Z-Rotations.

    """
    p = Program()

    if len(qubit_pairs) != gammas_pairs.declared_size:
        raise ValueError("zz_rotation_angles must have the same length as qubits_pairs")

    for gamma_pair, qubit_pair in zip(gammas_pairs, qubit_pairs):
        p.inst(RZ(2 * gamma_pair, qubit_pair[0]))
        p.inst(RZ(2 * gamma_pair, qubit_pair[1]))
        p.inst(CPHASE(-4 * gamma_pair, qubit_pair[0], qubit_pair[1]))

    if gammas_singles.declared_size != len(qubit_singles):
        raise ValueError("z_rotation_angles must have the same length as qubit_singles")

    for gamma_single, qubit in zip(gammas_singles, qubit_singles):
        p.inst(RZ(2 * gamma_single, qubit))

    return p

def smallish_ansatz(params):
    """Returns an ansatz with two parameters."""
    return Program(RX(params[0], 0), RZ(params[1], 0))

alpha_label = f"{prefix}_{suffix_alpha}"
    beta_label = f"{prefix}_{suffix_beta}"
    gamma_label = f"{prefix}_{suffix_gamma}"

    p = Program()

    alpha = p.declare(alpha_label, "REAL", len(qubits))
    beta = p.declare(beta_label, "REAL", len(qubits))
    gamma = p.declare(gamma_label, "REAL", len(qubits))

    for idx, q in enumerate(qubits):
        p += RZ(alpha[idx], q)
        p += RX(np.pi / 2, q)
        p += RZ(beta[idx], q)
        p += RX(-np.pi / 2, q)
        p += RZ(gamma[idx], q)

    return p