How to use ffsim’s Qiskit Sampler primitive

In this guide, we show how to use FfsimSampler, ffsim’s implementation of the Qiskit Sampler primitive, to sample from quantum circuits built from fermionic gates.

Example of using FfsimSampler

We’ll begin with a simple example of using the ffsim Sampler. First, let’s create an example circuit using gates from the ffsim.qiskit module.

[1]:
import numpy as np
from qiskit.circuit import QuantumCircuit, QuantumRegister

import ffsim

# Let's use 8 spatial orbitals with 4 alpha electrons and 4 beta electrons.
norb = 8
nelec = (4, 4)
n_alpha, n_beta = nelec

# Generate some random data
rng = np.random.default_rng(12345)
orbital_rotation = ffsim.random.random_unitary(norb, seed=rng)
diag_coulomb_mat = ffsim.random.random_real_symmetric_matrix(norb, seed=rng)

# Create an example circuit
qubits = QuantumRegister(2 * norb)
circuit = QuantumCircuit(qubits)
circuit.append(
    ffsim.qiskit.PrepareSlaterDeterminantJW(
        norb,
        occupied_orbitals=[range(n_alpha), range(n_beta)],
        orbital_rotation=orbital_rotation,
    ),
    qubits,
)
circuit.append(
    ffsim.qiskit.DiagCoulombEvolutionJW(norb, diag_coulomb_mat, time=1.0), qubits
)
circuit.append(ffsim.qiskit.OrbitalRotationJW(norb, orbital_rotation.T.conj()), qubits)
circuit.measure_all()

Now, we initialize the ffsim Sampler and use it to sample 10,000 shots from our circuit. The input to the Sampler is a list of primitive unified blocs, or PUBs. In the cell output we display only the top 10 most commonly encountered bitstrings.

[2]:
# Initialize ffsim Sampler
sampler = ffsim.qiskit.FfsimSampler(default_shots=10_000, seed=rng)

# Form PUB, submit job, retrieve job result, and extract first (and only) PUB result
pub = (circuit,)
job = sampler.run([pub])
result = job.result()
pub_result = result[0]

# Get counts
counts = pub_result.data.meas.get_counts()

# Display the 10 most commonly seen bitstrings and their counts
{k: v for k, v in sorted(counts.items(), key=lambda x: x[1], reverse=True)[:10]}
[2]:
{'0000111100001111': 210,
 '0000111100101011': 133,
 '0010101100001111': 101,
 '0000111100011101': 70,
 '0010011100101101': 66,
 '0010110100101011': 65,
 '0010110100100111': 57,
 '0001110100001111': 56,
 '0000111100011011': 48,
 '0010101100101101': 47}

Criteria for circuits that FfsimSampler can sample

FfsimSampler cannot sample from arbitrary Qiskit circuits. It can handle circuits that satisfy the following criteria:

Due to the need to preserve fermionic statistics, some of the supported gates must be applied to consecutive qubits in ascending order.

More examples

Sampling from an LUCJ circuit for a closed-shell molecule

The following code cell demonstrates a possible workflow for sampling from a spin-balanced LUCJ circuit for a nitrogen molecule in the 6-31g basis.

[3]:
import pyscf
import pyscf.cc
import pyscf.data.elements

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 4
active_space = range(n_frozen, mol.nao_nr())

# Get molecular data and Hamiltonian
scf = pyscf.scf.RHF(mol).run()
mol_data = ffsim.MolecularData.from_scf(scf, active_space=active_space)
norb, nelec = mol_data.norb, mol_data.nelec
mol_hamiltonian = mol_data.hamiltonian
print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()

# Use 2 ansatz layers
n_reps = 2
# Use interactions implementable on a square lattice
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = [(p, p) for p in range(norb)]
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    ccsd.t2, n_reps=n_reps, interaction_pairs=(pairs_aa, pairs_ab)
)

# Construct circuit
qubits = QuantumRegister(2 * norb)
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

# Sample 10,000 shots from the circuit using FfsimSampler
sampler = ffsim.qiskit.FfsimSampler(default_shots=10_000, seed=12345)
pub = (circuit,)
job = sampler.run([pub])
result = job.result()
pub_result = result[0]
counts = pub_result.data.meas.get_counts()

# Display the 10 most commonly seen bitstrings and their counts
{k: v for k, v in sorted(counts.items(), key=lambda x: x[1], reverse=True)[:10]}
converged SCF energy = -108.835236570775
norb = 14
nelec = (3, 3)
E(CCSD) = -108.9630419334855  E_corr = -0.1278053627110062
<class 'pyscf.cc.ccsd.CCSD'> does not have attributes  converged
[3]:
{'0000000000011100000000000111': 9924,
 '0000000000110100000000001101': 14,
 '0000000001110000000000000111': 10,
 '0000000000011100000000011100': 10,
 '0000000001011000000000010110': 9,
 '0001000001010000000000000111': 5,
 '0000000001011000100000000110': 4,
 '0100000000100100000000000111': 3,
 '0000000000011100100000001100': 3,
 '0010000000011000000000010110': 3}

Sampling from an LUCJ circuit for an open-shell molecule

The following code cell demonstrates a possible workflow for sampling from a spin-unbalanced LUCJ circuit for a hydroxyl radical in the 6-31g basis.

[4]:
import pyscf.data.elements
from pyscf import cc, gto

# Build HO molecule
mol = gto.Mole()
mol.build(
    atom=[["H", (0, 0, 0)], ["O", (0, 0, 1.1)]],
    basis="6-31g",
    spin=1,
    symmetry="Coov",
)

# Get molecular data and Hamiltonian
scf = pyscf.scf.ROHF(mol).run()
mol_data = ffsim.MolecularData.from_scf(scf)
norb, nelec = mol_data.norb, mol_data.nelec
mol_hamiltonian = mol_data.hamiltonian
print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = cc.CCSD(scf).run()

# Use 4 layers from opposite-spin amplitudes and 2 layers from same-spin amplitudes
n_reps = (4, 2)
# Use interactions implementable on a square lattice
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = [(p, p) for p in range(norb)]
pairs_bb = [(p, p + 1) for p in range(norb - 1)]
ucj_op = ffsim.UCJOpSpinUnbalanced.from_t_amplitudes(
    ccsd.t2, n_reps=n_reps, interaction_pairs=(pairs_aa, pairs_ab, pairs_bb)
)

# Construct circuit
qubits = QuantumRegister(2 * norb)
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
circuit.append(ffsim.qiskit.UCJOpSpinUnbalancedJW(ucj_op), qubits)
circuit.measure_all()

# Sample 10,000 shots from the circuit using FfsimSampler
sampler = ffsim.qiskit.FfsimSampler(default_shots=10_000, seed=12345)
pub = (circuit,)
job = sampler.run([pub])
result = job.result()
pub_result = result[0]
counts = pub_result.data.meas.get_counts()

# Display the 10 most commonly seen bitstrings and their counts
{k: v for k, v in sorted(counts.items(), key=lambda x: x[1], reverse=True)[:10]}
SCF not converged.
SCF energy = -75.3484557083881
norb = 11
nelec = (5, 4)

WARN: RCCSD method does not support ROHF method. ROHF object is converted to UHF object and UCCSD method is called.

E(UCCSD) = -75.45619739102422  E_corr = -0.1077416826361168
<class 'pyscf.cc.uccsd.UCCSD'> does not have attributes  converged
[4]:
{'0000000111100000011111': 9991,
 '0000100101100000111011': 2,
 '0000100110100000111011': 1,
 '0100000110100100001111': 1,
 '0101000001100000011111': 1,
 '0000010110100001011011': 1,
 '1000000101100000111011': 1,
 '0000000111100110000111': 1,
 '0000010101100001011011': 1}