Build Qiskit circuits for the LUCJ ansatz

This guide provides some examples of building and transpiling Qiskit circuits to implement the LUCJ ansatz.

import warnings

import pyscf
import pyscf.cc
from qiskit.circuit import QuantumCircuit, QuantumRegister
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import CouplingMap

import ffsim
from ffsim.qiskit import generate_lucj_pass_manager

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

LUCJ circuit for a closed-shell molecule

For a closed-shell system, use the spin-balanced LUCJ ansatz. This example creates an LUCJ circuit for a nitrogen molecule in the 6-31g basis and transpiled it for generic backend with square connectivity.

# 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 = 2
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
# Define interactions
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = (
    None  # Use generate_lucj_pass_manager() to find alpha-beta interactions later
)

# Create a generic backend with square connectivity
cmap_square = CouplingMap.from_grid(num_rows=12, num_columns=10)
backend = GenericBackendV2(
    num_qubits=cmap_square.size(),
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
    coupling_map=cmap_square,
)

# Generate a LUCJ-specific configured Qiskit preset pass manager
# and get pairs_ab, implementable on a square coupling map.
# The pass manager is confifgured with pre-init passes suggested by
# ffsim (ffsim.qiskit.PRE_INIT) and an error-aware LUCJ-friendly layout
pass_manager, pairs_ab = generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Use the backend implementable pairs_ab to construct the ucj_op
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    ccsd.t2,
    t1=ccsd.t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization.
    # Additionally, you may want to set the multi_stage_start or multi_stage_step
    # arguments (or both) to obtain a better result at increased computational cost.
    # See the API documentation for details.
    optimize=True,
    options=dict(maxiter=100),
)

# 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()

# Transpile the circuit with configured pass manager
transpiled = pass_manager.run(circuit)

transpiled.count_ops()

WARN: Unable to to identify input symmetry using original axes.
Different symmetry axes will be used.

converged SCF energy = -108.835236570774
norb = 16
nelec = (5, 5)
E(CCSD) = -109.0398256929732  E_corr = -0.2045891221988314

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12), (13, 13), (14, 14), (15, 15)].
Removing interaction (15, 15) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12), (13, 13), (14, 14)].
Removing interaction (14, 14) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12), (13, 13)].
Removing interaction (13, 13) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)].
Removing interaction (12, 12) from the end.

OrderedDict([('xx_plus_yy', 590),
             ('cp', 84),
             ('p', 64),
             ('measure', 32),
             ('x', 10),
             ('barrier', 1)])

As we initially set pairs_ab=None, the generate_lucj_pass_manager() starts with a default one. For square lattice, it is [(0, 0), (1, 1), (2, 2), ..., (norb-1, norb-1)]. However, if the hardware connectivity does not accommodate all starting pairs, then pairs are dropped from the end of the list until a valid layout is found, and a warning is raised. In the example above, last four pairs, from (15, 15) to (12, 12), are dropped with a warning message showing the dropped pairs.

Note that you can use a different pairs_ab during the ucj_op construction instead of the one returned by generate_lucj_pass_manager, but doing so may incur a higher SWAP overhead.

LUCJ circuit for an open-shell molecule

For an open-shell system, use the spin-unbalanced LUCJ ansatz. This example creates an LUCJ circuit for a hydroxyl radical in the 6-31g basis.

# Build HO molecule
mol = pyscf.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 = pyscf.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)]  # Update it later; can be ``None`` also
pairs_bb = [(p, p + 1) for p in range(norb - 1)]  # Can be ``None`` as well

# Generate the configured preset pass manager and implementable pairs_ab
pass_manager, pairs_ab = generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab, pairs_bb),
    optimization_level=3,
)

# Use the backend implementable pairs_ab to construct the ucj_op
ucj_op = ffsim.UCJOpSpinUnbalanced.from_t_amplitudes(
    ccsd.t2,
    t1=ccsd.t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab, pairs_bb),
    # Setting optimize=True enables the "compressed" factorization.
    # Additionally, you may want to set the multi_stage_start or multi_stage_step
    # arguments (or both) to obtain a better result at increased computational cost.
    # See the API documentation for details.
    optimize=True,
    options=dict(maxiter=100),
)

# 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()

# Transpile the circuit with configured pass manager
transpiled = pass_manager.run(circuit)

transpiled.count_ops()

SCF not converged.
SCF energy = -75.3484557085462
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.45619739099864  E_corr = -0.1077416824524348

OrderedDict([('xx_plus_yy', 718),
             ('p', 132),
             ('cp', 128),
             ('measure', 22),
             ('x', 9),
             ('barrier', 1)])

As norb=11 in this example, pairs_ab = [(0, 0), ..., (10, 10)], which the backend can accommodate. Therefore, this example does not raise warnings about dropped pairs.