python v1.x

Openfermion

Electronic structure package for quantum computers

Quick install

pip install openfermion openfermion-cirq

Background and History

OpenFermion was released in 2017 by a team at Google Quantum AI, with Ryan Babbush, Jarrod McClean, and Ian Kivlichan among the primary authors. The project addressed a specific gap in the quantum computing software stack: translating problems from quantum chemistry and materials science into qubit representations suitable for quantum algorithms. Before OpenFermion, researchers performing quantum chemistry calculations on quantum computers had to manually implement fermion-to-qubit mappings (Jordan-Wigner, Bravyi-Kitaev, and others), which was error-prone and repetitive. OpenFermion standardized these transformations into a reusable Python library.

The library’s core data types, FermionOperator and QubitOperator, provide symbolic algebra for fermionic and qubit Hamiltonians respectively. Users can construct molecular Hamiltonians from electronic structure calculations, apply qubit mappings, and generate the sparse matrices needed for classical verification or quantum algorithm input. The companion package openfermion-cirq provides circuit-level tools for implementing chemistry algorithms in Google’s Cirq framework, while openfermion-pyscf interfaces with the PySCF classical chemistry package for computing molecular integrals.

OpenFermion quickly became the standard open-source tool for quantum computational chemistry research. It was used in several of the early landmark quantum chemistry experiments on real hardware, including variational calculations of molecular ground states on Google’s processors. The library influenced the design of similar functionality in other frameworks, and both Qiskit Nature and PennyLane’s quantum chemistry modules address the same problem space that OpenFermion originally defined.

As of 2025, OpenFermion is maintained under the quantumlib GitHub organization alongside Cirq. It has over 1,500 GitHub stars and remains widely used in academic research on quantum algorithms for chemistry and materials science. Development activity has moderated compared to its early years, with the library considered mature and stable for its core functionality. Some researchers have migrated to Qiskit Nature or PennyLane’s chemistry tools for tighter integration with those ecosystems, but OpenFermion remains the reference implementation for fermion-to-qubit mappings and a common dependency in quantum chemistry research codebases.

Installation

pip install openfermion openfermion-cirq
pip install openfermion-pyscf   # for PySCF molecular integrals

Key Imports

from openfermion import FermionOperator, QubitOperator
from openfermion import jordan_wigner, bravyi_kitaev, get_sparse_operator
from openfermion.utils import hermitian_conjugated, count_qubits, normal_ordered

FermionOperator

^ denotes a creation operator; its absence denotes an annihilation operator.

from openfermion import FermionOperator
from openfermion.utils import hermitian_conjugated, normal_ordered

hop = FermionOperator('0^ 1', -0.5)           # a†_0 a_1 hopping
n_0 = FermionOperator('0^ 0',  1.0)           # number operator
sym = hop + hermitian_conjugated(hop)          # Hermitian hopping
print(normal_ordered(FermionOperator('1') * FermionOperator('0^')))
# -1.0 [0^ 1] + 1.0 []

Qubit Mappings

Jordan-Wigner (one qubit per spin orbital, O(n) Pauli weight):

from openfermion import jordan_wigner

op    = FermionOperator('0^ 1', 1.0) + FermionOperator('1^ 0', 1.0)
qb_op = jordan_wigner(op)
# -0.5j [X0 Y1] + 0.5j [Y0 X1]

Bravyi-Kitaev (binary-tree encoding, O(log n) Pauli weight, lower T-gate count):

from openfermion import bravyi_kitaev

qb_op_bk = bravyi_kitaev(op)

Molecule Workflow

from openfermion import MolecularData, jordan_wigner, get_sparse_operator
from openfermionpyscf import run_pyscf
import numpy as np

geometry  = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
molecule  = MolecularData(geometry, basis='sto-3g', multiplicity=1, charge=0)
molecule  = run_pyscf(molecule, run_fci=True)

qubit_h    = jordan_wigner(molecule.get_molecular_hamiltonian())
H_matrix   = get_sparse_operator(qubit_h).toarray()
fci_energy = np.linalg.eigvalsh(H_matrix)[0]
print(f"FCI energy: {fci_energy:.6f} Ha")  # -1.136190 Ha

QubitOperator and Sparse Matrix

from openfermion import QubitOperator, get_sparse_operator

H_q    = QubitOperator('Z0', 0.5) + QubitOperator('X0 X1', 0.25)
sparse = get_sparse_operator(H_q)      # scipy sparse
dense  = sparse.toarray()              # numpy array

InteractionOperator and Molecular Hamiltonians

InteractionOperator stores a molecular Hamiltonian in its most compact form: a constant (nuclear repulsion energy), a matrix of one-body integrals, and a tensor of two-body integrals. This is the natural output of a Hartree-Fock calculation and the standard input for fermion-to-qubit mappings.

from openfermion import InteractionOperator, get_fermion_operator, jordan_wigner
import numpy as np

# Nuclear repulsion energy
constant = 0.7151043

# One-body integrals (h_pq): kinetic energy + electron-nuclear attraction
# For H2 in STO-3G, this is a 4x4 matrix (4 spin orbitals)
one_body = np.array([
    [-1.2524, 0.0000, 0.0000, 0.0000],
    [ 0.0000, -1.2524, 0.0000, 0.0000],
    [ 0.0000,  0.0000, -0.4759, 0.0000],
    [ 0.0000,  0.0000,  0.0000, -0.4759],
])

# Two-body integrals (h_pqrs): electron-electron repulsion
# Stored in physicist's notation: <pq|rs>
two_body = np.zeros((4, 4, 4, 4))
two_body[0, 1, 1, 0] = 0.6745
two_body[1, 0, 0, 1] = 0.6745
two_body[2, 3, 3, 2] = 0.6636
two_body[3, 2, 2, 3] = 0.6636
two_body[0, 2, 2, 0] = 0.6631
two_body[2, 0, 0, 2] = 0.6631
two_body[1, 3, 3, 1] = 0.6631
two_body[3, 1, 1, 3] = 0.6631
two_body[0, 3, 3, 0] = 0.6973
two_body[3, 0, 0, 3] = 0.6973
two_body[1, 2, 2, 1] = 0.6973
two_body[2, 1, 1, 2] = 0.6973

hamiltonian = InteractionOperator(constant, one_body, two_body)

# Convert to FermionOperator for symbolic manipulation
fermion_h = get_fermion_operator(hamiltonian)
print(f"Fermion terms: {len(fermion_h.terms)}")

# Convert to qubit representation
qubit_h = jordan_wigner(fermion_h)
print(f"Qubit terms:   {len(qubit_h.terms)}")

In practice, you rarely build these arrays by hand. The MolecularData class populates them from a classical chemistry backend:

from openfermion import MolecularData
from openfermionpyscf import run_pyscf

geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
mol = MolecularData(geometry, basis='sto-3g', multiplicity=1, charge=0)
mol = run_pyscf(mol, run_fci=True)

# get_molecular_hamiltonian() returns an InteractionOperator
int_op = mol.get_molecular_hamiltonian()
print(f"Nuclear repulsion: {int_op.constant:.6f}")
print(f"One-body shape:    {int_op.one_body_tensor.shape}")
print(f"Two-body shape:    {int_op.two_body_tensor.shape}")

Active Space Approximation

Full orbital spaces grow quickly with basis set size. For a molecule with N spatial orbitals and M electrons, the qubit Hamiltonian requires 2N qubits (one per spin orbital). Active space methods reduce this by freezing core orbitals that are always occupied and discarding high-energy virtual orbitals that remain empty. Only the “active” frontier orbitals near the Fermi level are retained.

from openfermion import MolecularData, jordan_wigner, get_sparse_operator
from openfermionpyscf import run_pyscf
import numpy as np

# LiH in 6-31G: 11 spatial orbitals = 22 spin orbitals = 22 qubits
geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.6))]
mol = MolecularData(geometry, basis='6-31g', multiplicity=1, charge=0)
mol = run_pyscf(mol, run_fci=True)

full_hamiltonian = mol.get_molecular_hamiltonian()
n_qubits_full = 2 * mol.n_orbitals
print(f"Full space: {n_qubits_full} qubits ({mol.n_orbitals} spatial orbitals)")

# Reduce to an active space: 2 electrons in 2 spatial orbitals (4 qubits)
# occupied_indices: spatial orbitals frozen as doubly occupied
# active_indices:   spatial orbitals kept in the active space
occupied_indices = list(range(1))     # freeze orbital 0 (Li 1s core)
active_indices   = list(range(1, 3))  # keep orbitals 1 and 2

active_hamiltonian = mol.get_molecular_hamiltonian(
    occupied_indices=occupied_indices,
    active_indices=active_indices
)
qubit_h_active = jordan_wigner(active_hamiltonian)
n_qubits_active = 2 * len(active_indices)
print(f"Active space: {n_qubits_active} qubits ({len(active_indices)} spatial orbitals)")

# Verify: diagonalize the active-space Hamiltonian
H_sparse = get_sparse_operator(qubit_h_active)
eigenvalues = np.linalg.eigvalsh(H_sparse.toarray())
print(f"Active space ground energy: {eigenvalues[0]:.6f} Ha")
print(f"Full FCI energy:            {mol.fci_energy:.6f} Ha")

The occupied_indices and active_indices parameters on get_molecular_hamiltonian() handle the integral transformation internally. The frozen-core orbitals contribute a constant energy shift that is folded into the Hamiltonian’s constant term.

Trotterization

Hamiltonian simulation requires implementing the time-evolution operator exp(-iHt). Since Hamiltonian terms generally do not commute, the matrix exponential cannot be split into a product of individual term exponentials exactly. Trotter-Suzuki decomposition approximates this product with controllable error.

OpenFermion provides trotterize utilities through the openfermion-cirq companion package. You can also construct Trotter steps manually from the qubit Hamiltonian.

Manual Trotter step construction:

from openfermion import jordan_wigner, FermionOperator
from openfermion import get_sparse_operator
import numpy as np
from scipy.linalg import expm

# Simple two-site Hubbard model
t_hop = 1.0   # hopping strength
U     = 2.0   # on-site repulsion

hopping   = FermionOperator('0^ 1', -t_hop) + FermionOperator('1^ 0', -t_hop)
hopping  += FermionOperator('2^ 3', -t_hop) + FermionOperator('3^ 2', -t_hop)
repulsion = FermionOperator('0^ 0 2^ 2', U) + FermionOperator('1^ 1 3^ 3', U)

H_total = hopping + repulsion
H_qubit = jordan_wigner(H_total)
H_mat   = get_sparse_operator(H_qubit).toarray()

# Exact time evolution for comparison
time = 0.5
U_exact = expm(-1j * H_mat * time)

# First-order Trotter: exp(-iHt) ≈ exp(-iH_hop*t) * exp(-iH_rep*t)
H_hop_mat = get_sparse_operator(jordan_wigner(hopping)).toarray()
H_rep_mat = get_sparse_operator(jordan_wigner(repulsion)).toarray()

n_steps = 10
dt = time / n_steps
U_trotter = np.eye(H_mat.shape[0], dtype=complex)
for _ in range(n_steps):
    U_trotter = expm(-1j * H_hop_mat * dt) @ expm(-1j * H_rep_mat * dt) @ U_trotter

# Measure Trotter error via operator norm
error = np.linalg.norm(U_exact - U_trotter, ord=2)
print(f"Trotter error (n={n_steps}): {error:.6e}")

Second-order (symmetric) Trotter cuts the error from O(dt^2) to O(dt^3) per step:

# Second-order Trotter: exp(-iH_hop*dt/2) * exp(-iH_rep*dt) * exp(-iH_hop*dt/2)
U_trotter_2 = np.eye(H_mat.shape[0], dtype=complex)
for _ in range(n_steps):
    U_trotter_2 = (
        expm(-1j * H_hop_mat * dt / 2)
        @ expm(-1j * H_rep_mat * dt)
        @ expm(-1j * H_hop_mat * dt / 2)
        @ U_trotter_2
    )

error_2 = np.linalg.norm(U_exact - U_trotter_2, ord=2)
print(f"Trotter error (2nd order, n={n_steps}): {error_2:.6e}")

For circuit-level Trotterization that produces Cirq gates directly, see the openfermion-cirq package’s trotter module, which provides TrotterAlgorithm and TrotterStep classes.

Full VQE Workflow with Cirq

This end-to-end example builds the H2 Hamiltonian, maps it to qubits, constructs a hardware-efficient ansatz in Cirq, and runs a variational optimization loop on a Cirq simulator.

import numpy as np
import cirq
from scipy.optimize import minimize
from openfermion import MolecularData, jordan_wigner, get_sparse_operator
from openfermionpyscf import run_pyscf

# --- Step 1: Build the molecular Hamiltonian ---
geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
mol = MolecularData(geometry, basis='sto-3g', multiplicity=1, charge=0)
mol = run_pyscf(mol, run_fci=True)

qubit_hamiltonian = jordan_wigner(mol.get_molecular_hamiltonian())
H_sparse = get_sparse_operator(qubit_hamiltonian)
n_qubits = 4  # H2 in STO-3G: 2 spatial orbitals = 4 spin orbitals

print(f"Target FCI energy: {mol.fci_energy:.6f} Ha")
print(f"Hamiltonian terms: {len(qubit_hamiltonian.terms)}")

# --- Step 2: Define qubits and hardware-efficient ansatz ---
qubits = cirq.LineQubit.range(n_qubits)

def hardware_efficient_ansatz(params):
    """Single layer of Ry rotations followed by a CNOT entangling ladder."""
    circuit = cirq.Circuit()

    # Initial state: Hartree-Fock |0011> (2 electrons in lowest orbitals)
    circuit.append(cirq.X(qubits[0]))
    circuit.append(cirq.X(qubits[1]))

    # Layer of parameterized Ry rotations
    for i in range(n_qubits):
        circuit.append(cirq.ry(params[i]).on(qubits[i]))

    # CNOT entangling ladder
    for i in range(n_qubits - 1):
        circuit.append(cirq.CNOT(qubits[i], qubits[i + 1]))

    # Second layer of Ry rotations
    for i in range(n_qubits):
        circuit.append(cirq.ry(params[n_qubits + i]).on(qubits[i]))

    return circuit

n_params = 2 * n_qubits  # 8 parameters for this ansatz

# --- Step 3: Expectation value calculation ---
simulator = cirq.Simulator()

def expectation_value(params):
    """Compute <psi(params)|H|psi(params)> via statevector simulation."""
    circuit = hardware_efficient_ansatz(params)
    result = simulator.simulate(circuit)
    state = result.final_state_vector

    # Compute expectation using the sparse Hamiltonian matrix
    energy = np.real(state.conj() @ H_sparse @ state)
    return energy

# --- Step 4: Classical optimization loop ---
np.random.seed(42)
initial_params = np.random.uniform(-0.1, 0.1, n_params)

print(f"\nInitial energy:  {expectation_value(initial_params):.6f} Ha")

result = minimize(
    expectation_value,
    initial_params,
    method='COBYLA',
    options={'maxiter': 500, 'rhobeg': 0.5}
)

print(f"VQE energy:      {result.fun:.6f} Ha")
print(f"FCI energy:      {mol.fci_energy:.6f} Ha")
print(f"Error:           {abs(result.fun - mol.fci_energy):.6f} Ha")
print(f"Converged:       {result.success}")

This example uses statevector simulation for computing the exact expectation value. On real hardware or shot-based simulators, you would instead decompose the Hamiltonian into measurable Pauli strings and estimate each term’s expectation value from measurement statistics. The openfermion-cirq package provides utilities for this decomposition.

OpenFermion vs. Alternatives

OpenFermion, Qiskit Nature, and PennyLane’s quantum chemistry module all address the same core problem (fermion-to-qubit Hamiltonian construction for quantum algorithms), but they differ in design philosophy and integration:

  • Level of abstraction. OpenFermion exposes lower-level primitives (FermionOperator, InteractionOperator, individual mapping functions) and leaves circuit construction to the user or to openfermion-cirq. Qiskit Nature and PennyLane wrap more of the pipeline into high-level solver classes that handle mapping, ansatz construction, and optimization together.

  • Framework coupling. OpenFermion is framework-agnostic at its core, with Cirq integration provided by a separate package. Qiskit Nature is tightly coupled to Qiskit’s circuit and algorithm infrastructure. PennyLane’s chemistry tools are built into the PennyLane ecosystem and produce PennyLane-native Hamiltonians.

  • Autodifferentiation. PennyLane supports automatic differentiation of quantum circuits natively, which makes gradient-based VQE optimization straightforward. OpenFermion and Qiskit Nature rely on finite-difference gradients or parameter-shift rules implemented separately.

  • Maintenance and ecosystem. OpenFermion is mature and stable but receives less active development than in its early years. Qiskit Nature moved to community maintenance in 2024, with its long-term trajectory dependent on contributors. PennyLane’s chemistry module is actively developed by Xanadu and benefits from tight integration with their hardware and software roadmap.

  • Research vs. production. OpenFermion remains the most commonly cited library in quantum chemistry research papers and is well-suited to exploratory work where you want direct control over every transformation step. Qiskit Nature and PennyLane are better choices when you want a batteries-included workflow that handles circuit construction and execution end-to-end.

Learning Resources