Build an LUCJ ansatz

Important

The concepts in this guide are currently available only in the Python API. Equivalent functionality will be made available through the C API in a future release.

The local unitary cluster Jastrow (LUCJ) ansatz is a compact, hardware-efficient parametrization of a correlated electronic wavefunction. It is a member of the more general unitary cluster Jastrow (UCJ) family and takes the form

\[\lvert \Psi \rangle = \left(\prod_{k=1}^{L} \mathcal{U}_k\, e^{i \mathcal{J}_k}\, \mathcal{U}_k^\dagger\right) \lvert \Phi_0 \rangle,\]

where \(\lvert \Phi_0 \rangle\) is a reference state (typically Hartree-Fock), each \(\mathcal{U}_k\) is an orbital rotation, and each \(\mathcal{J}_k\) is a diagonal Coulomb operator

\[\mathcal{J} = \frac12 \sum_{ij,\sigma\tau} \mathbf{J}^{\sigma\tau}_{ij}\, n_{i\sigma}\, n_{j\tau},\]

with \(n_{i\sigma}\) the number operator on spatial orbital \(i\) with spin \(\sigma\). This guide shows how to assemble such an ansatz for a real molecule using the UCJ gate from qiskit_fermions.circuit.library.

1. Run the classical calculation

The (L)UCJ ansatz can be initialized from the amplitudes of a coupled-cluster singles and doubles (CCSD) calculation. Here we run restricted Hartree-Fock followed by CCSD for a hydrogen molecule in the 6-31g basis, using PySCF for the quantum chemistry.

>>> import pyscf
>>> import pyscf.cc
>>>
>>> # build the molecule and run Hartree-Fock
>>> mol = pyscf.gto.Mole()
>>> mol.build(
...     atom=[["H", (0, 0, 0)], ["H", (0, 0, 0.74)]],
...     basis="6-31g",
...     symmetry="Dooh",
...     verbose=0,
... )
<pyscf.gto.mole.Mole object at ...>
>>> scf = pyscf.scf.RHF(mol).run()
>>>
>>> mo_coeff = scf.mo_coeff
>>> norb = mo_coeff.shape[1]
>>> nelec = (mol.nelec[0], mol.nelec[1])
>>>
>>> # run CCSD for the t-amplitudes
>>> ccsd = pyscf.cc.CCSD(scf).run()
>>> t1, t2 = ccsd.t1, ccsd.t2

2. Build the molecular Hamiltonian as a fermionic operator

We will need the Hamiltonian later to evaluate the ansatz energy. We build it directly as a FermionOperator from the molecular-orbital integrals: the one-body integrals h1e (the core Hamiltonian in the MO basis), the two-body integrals h2e (from pyscf.ao2mo()), and the constant nuclear-repulsion energy. The electronic-integral constructors from_1body_tril_spin_sym() and from_2body_tril_spin_sym() expect the integrals in packed (lower-triangular) chemist ordering, which is exactly what PySCF produces.

>>> from pyscf import ao2mo, lib
>>>
>>> from qiskit_fermions.operators import FermionOperator
>>>
>>> # one- and two-body molecular-orbital integrals and the nuclear-repulsion energy
>>> h1e = mo_coeff.T @ scf.get_hcore() @ mo_coeff
>>> h2e = ao2mo.kernel(mol, mo_coeff)
>>> ecore = mol.energy_nuc()
>>>
>>> # pack into the lower-triangular chemist-ordered layout the constructors expect
>>> h1e_tril = lib.pack_tril(h1e)
>>> h2e_tril = lib.pack_tril(h2e)
>>>
>>> hamiltonian = ecore * FermionOperator.one()
>>> hamiltonian += FermionOperator.from_1body_tril_spin_sym(h1e_tril, norb)
>>> hamiltonian += FermionOperator.from_2body_tril_spin_sym(h2e_tril, norb)

// Given the packed integrals (obtained as above), the C API exposes matching
// electronic-integral constructors:
//
//   QfFermionOperator *h1 = qf_ferm_op_from_1body_tril_spin_sym(h1e_tril, norb);
//   QfFermionOperator *h2 = qf_ferm_op_from_2body_tril_spin_sym(h2e_tril, norb);

3. Build the LUCJ circuit

The UCJ gate assembles the ansatz directly from the coupled-cluster amplitudes. Its from_t_amplitudes() constructor performs a double factorization of the \(t_2\) amplitudes (via double_factorized_t2()) to obtain the per-layer diagonal Coulomb matrices and orbital rotations, and derives an optional final orbital rotation from the \(t_1\) amplitudes.

The number of ansatz repetitions \(L\) equals the number of terms in the double factorization. Truncating it with the n_reps argument trades some accuracy for a shallower circuit; here we keep the two largest terms, which recovers most of the correlation energy while halving the number of layers.

>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import UCJ
>>>
>>> ansatz = UCJ.from_t_amplitudes(nelec, t2, t1=t1, n_reps=2)
>>>
>>> circuit = FermionicCircuit(2 * norb)
>>> circuit.append(ansatz, circuit.modes)

// The C API for FermionicCircuit is not available yet.

Decomposing the gate reveals its anatomy: an InitializeModes gate prepares the Hartree-Fock reference determinant, and each ansatz layer contributes an OrbitalRotation \(\mathcal{U}_k^\dagger\), then \(e^{i\mathcal{J}_k}\) (an Evolution of the diagonal Coulomb operator \(\mathcal{J}_k\)), then \(\mathcal{U}_k\), with a final OrbitalRotation at the end. The orbital rotations act per spin sector, so each is placed on the alpha modes 0..norb and the beta modes norb..2*norb independently.

>>> circuit.decompose().draw("mpl", fold=-1)
<Figure size ... with 1 Axes>

(png, hires.png, pdf)

The gates that the UCJ ansatz decomposes into.

Note

Each layer ends with \(\mathcal{U}_k\) and the next begins with \(\mathcal{U}_{k+1}^\dagger\), so adjacent OrbitalRotation gates could be merged into a single rotation. A transpilation pass performing this fusion is a planned future development.

4. Simulate the ansatz and evaluate its energy

Because every gate in the circuit implements ffsim’s ffsim.SupportsApplyUnitary protocol, the whole FermionicCircuit can be applied to a fixed particle-number state vector with ffsim.apply_unitary(), starting from the Hartree-Fock reference. The FermionOperator likewise implements ffsim’s ffsim.SupportsLinearOperator protocol, so we can obtain a SciPy LinearOperator for it via ffsim.linear_operator() and evaluate the ansatz energy as the expectation value of the molecular Hamiltonian.

>>> import ffsim
>>> import numpy as np
>>>
>>> reference = ffsim.hartree_fock_state(norb, nelec)
>>> state = ffsim.apply_unitary(reference, circuit, norb=norb, nelec=nelec)
>>>
>>> linop = ffsim.linear_operator(hamiltonian, norb=norb, nelec=nelec)
>>> energy = np.vdot(state, linop @ state).real
>>> print(f"LUCJ energy: {energy:.8f} Hartree")
LUCJ energy: -1.14618323 Hartree

// Circuit simulation is not available via the C API.

The LUCJ energy improves substantially on the Hartree-Fock reference and approaches the CCSD energy it was initialized from – the small remaining gap is the price of truncating the ansatz to two repetitions:

>>> print(f"Hartree-Fock: {scf.e_tot:.8f} Hartree")
Hartree-Fock: -1.12675532 Hartree
>>> print(f"CCSD:         {ccsd.e_tot:.8f} Hartree")
CCSD:         -1.15167268 Hartree

Note

The UCJ gate is initialized here from an exact double factorization of the \(t_2\) amplitudes. ffsim additionally offers an optimized (“compressed”) double factorization (its from_t_amplitudes(..., optimize=True)), which has no equivalent in this package. To use it, build an ffsim UCJ operator with optimize=True and pass its diag_coulomb_mats / orbital_rotations / final_orbital_rotation into UCJ directly.

Next steps