Building Quantum Programs with PyQuil and Quilbase

PyQuil is Rigetti’s Python SDK for writing quantum programs in Quil (Quantum Instruction Language). Starting with version 4.x, PyQuil reorganized its API around the quilbase module, which provides the low-level building blocks for programs: instructions, classical memory declarations, and gate objects. This tutorial walks through constructing, parameterizing, and executing programs using the modern API. Along the way, you will learn how to compose modular circuits, define custom gates, sweep parameters for variational algorithms, and inspect what the compiler produces.

Installation and Setup

pip install pyquil

PyQuil 4.x execution against real Rigetti hardware requires a QCS account and credentials. For local simulation, use the QVM (quantum virtual machine) and quilc compiler, both available as Docker images:

docker run --rm -it -p 5000:5000 rigetti/qvm -S
docker run --rm -it -p 6000:6000 rigetti/quilc -R

With both services running, PyQuil can compile and simulate programs locally without any QCS account. The QVM is a full wavefunction simulator, so it produces exact results (up to shot noise from sampling).

What Quilbase Is and Why It Exists

When you write a quantum program, one approach is to construct Quil source code as a raw string:

quil_string = """
DECLARE ro BIT [2]
H 0
CNOT 0 1
MEASURE 0 ro[0]
MEASURE 1 ro[1]
"""

This works, but it is error-prone and hard to manipulate programmatically. A typo like CONT 0 1 instead of CNOT 0 1 only surfaces when the compiler parses the string. You also cannot easily iterate over the instructions, filter for specific gates, or compose two programs together.

Quilbase solves these problems by modeling every Quil construct as a Python object. The Program class holds a list of these objects rather than a flat string. This design gives you three concrete benefits:

Type checking. If you misspell a gate name, you get an immediate Python error rather than a cryptic compiler message later:

from pyquil.quilbase import Gate, Qubit

# This works: "H" is a valid gate name
h_gate = Gate("H", [], [Qubit(0)])

# This also "works" at construction time, but quilc will reject it.
# The advantage is you can validate gate names in Python before sending to quilc.

Introspection. You can iterate over a program’s instructions and inspect them:

from pyquil import Program
from pyquil.gates import H, CNOT, MEASURE
from pyquil.quilbase import Gate, Measurement

p = Program(H(0), CNOT(0, 1))
ro = p.declare("ro", "BIT", 2)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

for instr in p.instructions:
    if isinstance(instr, Gate):
        print(f"Gate: {instr.name}, qubits: {instr.qubits}")
    elif isinstance(instr, Measurement):
        print(f"Measure qubit {instr.qubit}")

Composition. You can merge programs using the + operator, which concatenates their instruction lists:

from pyquil import Program
from pyquil.gates import H, CNOT

a = Program(H(0))
b = Program(CNOT(0, 1))
combined = a + b
print(combined)
# H 0
# CNOT 0 1

These properties make Quilbase the right foundation for any serious PyQuil work, from hand-crafted experiments to auto-generated variational circuits.

The Program Object

The central object in PyQuil is Program. It accumulates instructions and tracks declared classical memory regions.

from pyquil import Program
from pyquil.gates import H, CNOT, MEASURE
from pyquil.quilbase import Declare

# Declare a classical memory register
p = Program()
ro = p.declare("ro", "BIT", 2)

# Build a Bell state
p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

print(p)

The Declare instruction is important: in Quil, classical memory must be declared before use. The declare helper method on Program creates the declaration and returns a MemoryReference you can pass directly to MEASURE.

You can also inspect useful properties of a program:

# Get all qubit indices used in the program
print(p.get_qubits())  # [0, 1]

# Get all memory declarations
print(p.declarations)  # {'ro': Declare('ro', 'BIT', 2)}

Program Composition

One of the most powerful patterns in Quilbase is building modular circuit components and composing them. You define reusable building blocks, then combine them into a full program:

from pyquil import Program
from pyquil.gates import H, CNOT, MEASURE

# State preparation block
state_prep = Program()
state_prep += H(0)
state_prep += CNOT(0, 1)

# Measurement block
measurement = Program()
ro = measurement.declare("ro", "BIT", 2)
measurement += MEASURE(0, ro[0])
measurement += MEASURE(1, ro[1])

# Compose the full program
full_program = state_prep + measurement
print(full_program)

This prints:

H 0
CNOT 0 1
DECLARE ro BIT[2]
MEASURE 0 ro[0]
MEASURE 1 ro[1]

Composition is especially useful in variational algorithms where the ansatz and measurement routines are defined separately and swapped independently. You can also use += to append instructions or entire programs in place:

p = Program()
p += state_prep
p += measurement

Classical Memory and MemoryReference

Quil supports four classical memory types. Each has a specific size and purpose:

• BIT: 1-bit value (0 or 1). Used for measurement results. This is the type you use most often.
• OCTET: 8-bit unsigned integer. Used for byte-level data. Rarely needed in typical quantum programs.
• INTEGER: 64-bit signed integer. Used for loop counters or classical control flow indices in more advanced Quil programs.
• REAL: 64-bit floating point. Used for rotation angles in parameterized circuits. This is the key type for VQE and QAOA parameter sweeps, because you can patch REAL memory at runtime without recompiling.

Here is a program that declares all four types:

from pyquil import Program
from pyquil.gates import RZ, RX, MEASURE

p = Program()

# Measurement results
ro = p.declare("ro", "BIT", 2)

# Rotation angles for parameterized gates
angles = p.declare("angles", "REAL", 2)

# A counter for classical control flow
counter = p.declare("counter", "INTEGER", 1)

# A byte register (rarely used, shown for completeness)
data = p.declare("data", "OCTET", 1)

# Use REAL memory as gate parameters
p += RZ(angles[0], 0)
p += RX(angles[1], 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

print(p)

The angles[0] reference is a MemoryReference object, not a Python float. When you later execute the program, you supply the actual values through a memory map. This distinction is critical: using a MemoryReference lets the compiler produce a single compiled executable that you can run many times with different parameter values. Using a Python float bakes the value into the program, requiring recompilation for each new value.

Parameterized Execution with QVM

To run parameterized programs on the local QVM, use the get_qc interface:

from pyquil import Program, get_qc
from pyquil.gates import H, RZ, MEASURE

p = Program()
ro = p.declare("ro", "BIT", 1)
theta = p.declare("theta", "REAL", 1)

p += H(0)
p += RZ(theta[0], 0)
p += MEASURE(0, ro[0])
p.wrap_in_numshots_loop(1000)

# Get a virtual quantum computer (QVM-backed)
qc = get_qc("1q-qvm")

# Execute with a specific value for theta
results = qc.run(
    qc.compile(p),
    memory_map={"theta": [0.5]}   # radians
).readout_data["ro"]

ones = sum(r[0] for r in results)
print(f"P(|1>) estimate: {ones / 1000:.3f}")

The wrap_in_numshots_loop(1000) call is essential. Without it, the QVM runs only a single shot and you get a single measurement outcome instead of statistics.

VQE-Style Parameter Sweep

A common pattern in variational quantum algorithms is to compile a parameterized circuit once, then sweep over many parameter values. This avoids calling quilc repeatedly, which is critical for performance: each qc.compile() invocation runs the full compilation pipeline (gate decomposition, qubit routing, optimization passes). Compiling once and patching parameters at runtime can be 20x faster than recompiling for each value.

Here is a practical example that sweeps a rotation angle and computes the expectation value of the ZZ operator:

from pyquil import Program, get_qc
from pyquil.gates import RY, CNOT, MEASURE
import numpy as np

# Build a parameterized ansatz
p = Program()
theta = p.declare("theta", "REAL", 1)
ro = p.declare("ro", "BIT", 2)

p += RY(theta[0], 0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])
p.wrap_in_numshots_loop(500)

qc = get_qc("2q-qvm")

# Compile once
compiled = qc.compile(p)

# Sweep theta from 0 to 2*pi
energies = []
for theta_val in np.linspace(0, 2 * np.pi, 20):
    result = qc.run(compiled, memory_map={"theta": [theta_val]})
    counts = result.readout_data["ro"]
    # Compute expectation value of ZZ:
    # ZZ eigenvalue is +1 when both qubits agree, -1 when they disagree
    zz_vals = [(-1)**counts[i][0] * (-1)**counts[i][1] for i in range(500)]
    energies.append(np.mean(zz_vals))

# Print the sweep results
for angle, energy in zip(np.linspace(0, 2 * np.pi, 20), energies):
    print(f"theta={angle:.2f}  <ZZ>={energy:.3f}")

In a real VQE loop, you would replace the linear sweep with a classical optimizer (such as scipy.optimize.minimize) that chooses the next parameter values based on the measured energy.

Building Programs with Quilbase Directly

For more control, use quilbase objects directly to construct instructions without the gate helper functions.

from pyquil import Program
from pyquil.quilbase import (
    Gate,
    Declare,
    Measurement,
    MemoryReference,
    Qubit,
)
import numpy as np

p = Program()
theta = p.declare("theta", "REAL", 1)
ro = p.declare("ro", "BIT", 2)

# Construct a gate directly from quilbase
h_gate = Gate("H", [], [Qubit(0)])
cnot_gate = Gate("CNOT", [], [Qubit(0), Qubit(1)])

p += h_gate
p += cnot_gate
p += Measurement(Qubit(0), MemoryReference("ro", 0))
p += Measurement(Qubit(1), MemoryReference("ro", 1))

print(p)

This approach is useful when you generate programs programmatically, for example when building variational circuits in a loop or translating from another IR.

Inspecting Program Contents

Quilbase makes it straightforward to iterate over and inspect a program’s instructions. This is valuable for debugging, for writing custom transpilation passes, or for extracting circuit statistics:

from pyquil import Program
from pyquil.gates import H, CNOT, RZ, MEASURE
from pyquil.quilbase import Declare, Gate, Measurement

p = Program()
ro = p.declare("ro", "BIT", 2)
theta = p.declare("theta", "REAL", 1)

p += H(0)
p += RZ(theta[0], 0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

# Iterate over all instructions
for instr in p.instructions:
    if isinstance(instr, Gate):
        print(f"Gate: {instr.name}, qubits: {instr.qubits}")
    elif isinstance(instr, Measurement):
        print(f"Measure qubit {instr.qubit}")
    elif isinstance(instr, Declare):
        print(f"Declare: {instr.name} {instr.memory_type}[{instr.memory_size}]")

# Useful Program properties
print(f"Qubits used: {p.get_qubits()}")
print(f"Declarations: {p.declarations}")

Custom Gate Definitions with DefGate

Quilbase supports defining custom unitary gates using DefGate. This appends a DEFGATE block to the Quil program, making the gate available by name throughout the rest of the program.

import numpy as np
from pyquil import Program
from pyquil.gates import MEASURE
from pyquil.quilbase import DefGate

# Define a custom sqrt-X gate
sqrt_x_matrix = np.array([
    [0.5 + 0.5j, 0.5 - 0.5j],
    [0.5 - 0.5j, 0.5 + 0.5j]
])
sqrt_x_def = DefGate("SQRT-X", sqrt_x_matrix)

# get_constructor() returns a callable that produces gate instructions
SQRT_X = sqrt_x_def.get_constructor()

p = Program()
ro = p.declare("ro", "BIT", 1)

# Add the gate definition to the program (required before using the gate)
p += sqrt_x_def

# Apply the custom gate to qubit 0
p += SQRT_X(0)
p += MEASURE(0, ro[0])

print(p)

The printed output includes the full DEFGATE matrix definition followed by the gate application. When this program is compiled, quilc decomposes the custom gate into native hardware gates.

You can define multi-qubit custom gates the same way by providing a larger matrix. For a two-qubit gate, supply a 4x4 unitary matrix and call the constructor with two qubit arguments:

import numpy as np
from pyquil import Program
from pyquil.quilbase import DefGate

# Define a custom two-qubit gate (iSWAP as an example)
iswap_matrix = np.array([
    [1, 0, 0, 0],
    [0, 0, 1j, 0],
    [0, 1j, 0, 0],
    [0, 0, 0, 1]
])
iswap_def = DefGate("MY-ISWAP", iswap_matrix)
MY_ISWAP = iswap_def.get_constructor()

p = Program()
p += iswap_def
p += MY_ISWAP(0, 1)
print(p)

Parametric Gate Definitions with DefGateByPaulis

For variational quantum chemistry and VQE circuits, gates are often specified as exponentials of Pauli operators. DefGateByPaulis lets you define such gates directly without computing the matrix yourself:

from pyquil import Program
from pyquil.quilbase import DefGateByPaulis
from pyquil.quilatom import Parameter
from pyquil.paulis import PauliSum, PauliTerm

# Define a parameterized ZZ rotation gate
theta = Parameter("theta")
zz_gate_def = DefGateByPaulis(
    name="ZZ_ROTATION",
    parameters=[theta],
    paulis=PauliSum([PauliTerm("Z", 0) * PauliTerm("Z", 1)])
)

p = Program()
p += zz_gate_def
print(p)

This is the standard form for expressing time-evolution operators in quantum simulation. The compiler handles the decomposition into native gates. Using DefGateByPaulis rather than manually computing matrix exponentials reduces the chance of errors in your gate definitions and makes the program’s intent clearer.

Pragma Statements

Quilbase supports Quil PRAGMA directives, which communicate metadata to the compiler and simulator. Pragmas do not affect the quantum state directly; instead, they control how quilc compiles the program.

from pyquil import Program
from pyquil.gates import H, CNOT
from pyquil.quilbase import Pragma

p = Program()

# Prevent the compiler from reordering gates in this region
p += Pragma("PRESERVE_BLOCK")
p += H(0)
p += CNOT(0, 1)
p += Pragma("END_PRESERVE_BLOCK")

# Control how logical qubits map to physical qubits
p += Pragma("INITIAL_REWIRING", freeform_string='"PARTIAL"')

print(p)

The most commonly used pragmas are:

• PRESERVE_BLOCK / END_PRESERVE_BLOCK: Tells quilc not to reorder, cancel, or optimize gates within this region. Use this when you need exact control over gate ordering, for example in benchmarking or tomography circuits where reordering would change what you are measuring.

• INITIAL_REWIRING: Controls how logical qubits in your program map to physical qubits on the chip. The value "PARTIAL" lets the compiler choose an initial mapping but does not insert SWAP gates. The value "NAIVE" maps logical qubit N to physical qubit N with no remapping.

Connecting to Rigetti QCS

For execution on real Rigetti hardware, use the QCSClient and the get_qc function with a real device name:

from pyquil import Program, get_qc
from pyquil.gates import H, CNOT, MEASURE
from pyquil.api import QCSClient

# QCSClient reads credentials from ~/.qcs/credentials by default
client = QCSClient.load()

p = Program()
ro = p.declare("ro", "BIT", 2)
p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])
p.wrap_in_numshots_loop(100)

# Replace "Ankaa-2" with the device available on your QCS account
qc = get_qc("Ankaa-2", client=client)
compiled = qc.compile(p)
result = qc.run(compiled)
print(result.readout_data["ro"])

The compile step runs the Quil program through quilc, Rigetti’s open-source optimizing compiler, which rewrites the circuit into the hardware’s native gate set and maps logical qubits to physical qubits on the chip topology.

QCS Credentials and Authentication

The QCSClient looks for credentials in several places, checked in order:

1. The file ~/.qcs/credentials, created by running qcs configure from the command line.
2. Environment variables: QCS_SETTINGS_APPLICATIONS_PYQUIL_QVM_URL and QCS_SETTINGS_APPLICATIONS_PYQUIL_QUILC_URL control the QVM and quilc endpoints.
3. Programmatic configuration: you can pass auth_server and credentials arguments directly to QCSClient.

To verify that your credentials are working and see which QPUs are available:

from pyquil.api import QCSClient, list_quantum_computers

client = QCSClient.load()

# List available QPUs (not QVMs)
qpus = list_quantum_computers(client=client, qpus=True, qvms=False)
print([qpu.id for qpu in qpus])

Inspecting the Compiled Program

You can inspect the native Quil produced by quilc to understand what the compiler did:

from pyquil import Program, get_qc
from pyquil.gates import H, CNOT, MEASURE

qc = get_qc("2q-qvm")
p = Program()
ro = p.declare("ro", "BIT", 2)

p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])

compiled = qc.compile(p)
print(compiled.program)   # prints the native Quil after optimization

Understanding the Compiled Output

The compiled output contains several types of instructions that differ from your original program:

PRAGMA INITIAL_REWIRING “PARTIAL”: This tells you which qubit mapping strategy quilc chose. With "PARTIAL", the compiler found an initial assignment of logical qubits to physical qubits but did not insert additional SWAP operations. If you see "GREEDY", the compiler used a more aggressive strategy that may insert SWAPs to route two-qubit gates across the chip.

Native gate instructions: For the Rigetti Ankaa-2 processor, the native gate set consists of RZ, RX, and CZ. These correspond directly to the microwave and flux pulse operations the hardware can execute. Your H and CNOT gates are decomposed into sequences of these native gates. For example, an H gate typically becomes RZ(pi/2) RX(pi/2) RZ(pi/2).

PRAGMA PRESERVE_BLOCK / END_PRESERVE_BLOCK: These appear when the compiler marks regions it chose not to reorder. The gates inside these blocks execute in exactly the listed order.

Qubit index changes: Logical qubit 0 in your program may map to physical qubit 17 on the hardware. The compiled program uses physical qubit indices, and the pragmas at the top document the mapping. Always check these indices if you need to interpret raw results in terms of your original logical qubits.

Noise Simulation with QVM

Before committing to QPU time, you can estimate how noise will affect your circuit by using a noisy QVM:

from pyquil import Program, get_qc
from pyquil.gates import H, CNOT, MEASURE

# Get a noisy QVM that models Ankaa-2 noise characteristics
qc_noisy = get_qc("Ankaa-2", as_qvm=True, noisy=True)

p = Program()
ro = p.declare("ro", "BIT", 2)
p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])
p.wrap_in_numshots_loop(1000)

compiled = qc_noisy.compile(p)
result = qc_noisy.run(compiled)
counts = result.readout_data["ro"]

# Count Bell state outcomes
from collections import Counter
outcomes = Counter(tuple(row) for row in counts)
print(outcomes)

When you pass as_qvm=True, PyQuil creates a QVM simulator instead of connecting to the real QPU. Adding noisy=True loads a noise model calibrated to the real hardware’s gate error rates, T1 relaxation times, and T2 dephasing times. The results you get from this noisy simulation approximate what you would see on the actual device, including reduced fidelity and asymmetric error rates across different qubits.

Comparing noisy QVM results to noiseless QVM results (the default, with noisy=False) gives you a quick estimate of how much noise will impact your algorithm before you spend QPU credits.

Common Mistakes

Working with PyQuil and Quilbase, several mistakes come up repeatedly. Here is what to watch for:

Forgetting wrap_in_numshots_loop

Without wrap_in_numshots_loop(), the QVM runs exactly one shot. You get a single measurement outcome instead of the statistical distribution you almost certainly want:

# Wrong: only 1 shot
p = Program()
ro = p.declare("ro", "BIT", 1)
p += H(0)
p += MEASURE(0, ro[0])

qc = get_qc("1q-qvm")
result = qc.run(qc.compile(p))
print(result.readout_data["ro"])  # shape: (1, 1) -- just one sample

# Right: 1000 shots
p.wrap_in_numshots_loop(1000)
result = qc.run(qc.compile(p))
print(result.readout_data["ro"])  # shape: (1000, 1)

Using Python Floats Instead of MemoryReference

When you write RZ(0.5, 0), the float 0.5 gets baked into the program. This works but prevents you from patching the parameter at runtime. If you plan to sweep over values, always use a MemoryReference:

# Baked-in value: works but requires recompilation for each new angle
p = Program()
p += RZ(0.5, 0)

# Parameterized: compile once, run many times with different values
p = Program()
theta = p.declare("theta", "REAL", 1)
p += RZ(theta[0], 0)

Running a Program Without Compiling

You cannot pass a Program directly to qc.run(). You must compile first:

from pyquil import Program, get_qc
from pyquil.gates import H, MEASURE

p = Program()
ro = p.declare("ro", "BIT", 1)
p += H(0)
p += MEASURE(0, ro[0])
p.wrap_in_numshots_loop(100)

qc = get_qc("1q-qvm")

# Wrong: passing Program directly
# result = qc.run(p)  # TypeError

# Right: compile first, then run
compiled = qc.compile(p)
result = qc.run(compiled)

Misunderstanding Qubit Indices

Logical qubit 0 in your program does not necessarily correspond to physical qubit 0 on the hardware. The compiler maps logical qubits to physical qubits based on the chip topology and gate connectivity. On Ankaa-2, your logical qubit 0 might become physical qubit 17. The compiled program’s pragmas document this mapping. If you need specific physical qubits, use PRAGMA INITIAL_REWIRING "NAIVE" to force a direct mapping, though this may result in more SWAP operations.

Declaring Memory After Using It

Classical memory declarations must appear before any instructions that reference them. PyQuil generally handles ordering correctly when you use p.declare(), but if you construct instructions manually and add them out of order, you can end up with invalid Quil:

from pyquil import Program
from pyquil.gates import H, MEASURE

# Safe: declare first, then use
p = Program()
ro = p.declare("ro", "BIT", 1)
p += H(0)
p += MEASURE(0, ro[0])

Next Steps

From here, the natural extensions are writing full VQE or QAOA optimization loops that patch REAL memory at each optimizer step, defining problem-specific ansatze with DefGate and DefGateByPaulis, and using noisy QVM simulations to benchmark circuits before running on hardware. The PyQuil documentation covers the full API, and the Quil specification details every instruction the language supports.