Digital quantum simulation for spin and fermionic systems

The Hamiltonians classes

SpinHamiltonian, FermionHamiltonian and ElectronicStructureHamiltonian are essentially Observable classes with additional functionalities. Let us detail the additional features they contain.

Atomic and molecular studies study

Common many-body Hamiltonians

We present here Hamiltonian constructors for common many-body Hamiltonians.

The spin-fermion transforms

Fermionic ansatz circuits

Quantum phase estimation

qat.fermion.phase_estimation.perform_phase_estimation(H_el: ElectronicStructureHamiltonian, n_phase_bits: int, n_trotter_steps: int, init_vec: Optional[str] = None, n_adiab_steps: Optional[int] = 0, E_target: Optional[float] = 0, size_interval: Optional[float] = 2.0, basis_transform: Optional[str] = 'jordan-wigner', qpu=None, n_shots: Optional[int] = 0, verbose: Optional[bool] = False) → Tuple[float, float]

Perform quantum phase estimation (QPE) on an ElectronicStructureHamiltonian. This Hamiltonian is transformed to the computational basis via a Jordan-Wigner transformation and approximated via first order trotterization. Other transformations like parity and Bravyi-Kitaev are also possible.

When providing an initial state one can specify it either as a string composed of zeros and ones, or as a QRoutine which will produce it. The QPE is meant to start from an eigenstate of the Hamiltonian, however, knowing apriori even one eigenstate of the system may be challenging. Therefore, this function comes with adiabatic state preparation - an optional preliminary step to create an eigenstate of the Hamiltonian \(H\). This step consists in performing QPE n_adiab_steps number of times, but not to read the phase bits (it is set to only one), but rather to collapse the system to an eigenstate (read from the data bits). The first of the series of QPE executions starts from the lowest energy eigenstate of the Hamiltonian composed of \(h_{pq}\). Then, \(h_{pq}\) is linearly transformed to \(H\) and at each new step we start from the eigenstate of the Hamiltonian of the previous step. This guarantees that when the main QPE routine starts, it will do so from an eigenstate of the full \(H\).

Usually, energies lie outside the range \((-\frac{2\pi}{t}, 0)\). However, this range can be adjusted by searching inside the window \((E_{target} - \frac{\Delta}{2}, E_{target} + \frac{\Delta}{2})\) with \(E_{target}\) and \(\Delta\) specified by E_target and size_interval, respectively. It is suggested to always start from a large size interval and unbiased target energy like \(0\) thus enclosing many of the eigenenergies including the desired one. One can then narrow the window around an already found eigenenergy for a better precision. Working with a window not enclosing an eigenenergy would still evaluate to a result, but it may be misleading.

Warning

• Regarding the adiabatic state preparation, if the lowest energy eigenstate of the first-step Hamiltonian \(h_{pq}\) is also an eigenstate of the whole \(H\), the system will remain in it until the end of the whole adiabatic stage. Hence, this eigenstate may not be the one of the lowest energy anymore.

• As a rule of thumb, if small changes to the interval cause considerable deviations in the energy, that’s a sign that the window is too small or a different target energy may be better.

Parameters

• H_el (ElectronicStructureHamiltonian) – An electronic-structure Hamiltonian.

• n_phase_bits (int) – Number of qubits for the phase evaluation. The larger it is, the more accurate is the result.

• n_trotter_steps (int) – Number of first order trotterization steps. For good phase estimation it should also increase if n_phase_bits is increased.

• init_vec (Optional[str]) – Initial vector specified in the computational basis as a string - ‘01101’ for example. Starting from |0..0> an X will be applied to the respective qubits so as to produce the provided vector. This vector will enter the adiabatic state preparation routine if n_adiab_steps is not 0 or will be given straight to the main QPE routine.

• n_adiab_steps (Optional[int]) – Number of steps to pass from the part of the Hamiltonian containing only c_p^dagger * c_p terms (which is diagonal and fast to deal with) to the Hamiltonian of interest.

• E_target (Optional[float]) – Expected energy. If unknown, set to 0.

• size_interval (Optional[float]) – Size \(\Delta\) of the interval one thinks the value of the energy is in: \(E \in [E_\mathrm{target}-\Delta/2, E_\mathrm{target}+\Delta/2]\) If no idea take \(\Delta =2 E_\mathrm{max}\), with \(E_\mathrm{max}\) an upper bound of the energy.

• basis_transform (Optional[str]) – Transformation to go from qat.fermion.hamiltonians.ElectronicStructureHamiltonian into a qat.fermion.hamiltonians.SpinHamiltonian: one can use the “jordan-wigner” (default), “bravyi-kitaev” or “parity” transformations.

• qpu (Optional[QPU]) – QPU to use for computation. Will use by default the default installed QPU.

Returns

• Energy found,

• associated probability.

Return type

Tuple[float, float]

Note

Usually, energies lie outside the range \((-\frac{2\pi}{t}, 0)\). However, this range can be adjusted by specifying the arguments E_target and size_interval thus searching inside the window \((E_{t} - \frac{\Delta}{2}, E_{target} + \frac{size\_interval}{2})\), where \(E_{t}\) and \(\Delta\) stand for . We suggest to always start from a large size interval and unbiased target energy like 0 thus enclosing many of the eigenenergies including the desired one. One can then narrow the window around an already found eigenenergy for a better precision. Experience shows that working with a window not enclosing an eigenenergy makes the QPE still output a result, but it is misleading.

Quantum subspace expansion

qat.fermion.chemistry.qse.apply_quantum_subspace_expansion(hamiltonian: SpinHamiltonian, state_prep_circ: Circuit, expansion_operators: List[Observable], qpu: QPUHandler, nbshots: int = 0, threshold: float = 1e-15, return_matrices: bool = False) → Tuple[float, Optional[ndarray], Optional[ndarray]]

Apply quantum subspace expansion (QSE) to the given Hamiltonian.

QSE is a method that improves the quality of noisy results at the cost of additional measurements that help write the Hamiltonian in the small subspace where it can be diagonalized classically.

If \(\langle \Psi^\star | \hat{H} | \Psi^\star \rangle\) is the VQE result, the projected Hamiltonian matrix is built from

\[H^{(s)}_{i,j} = \langle \Psi^\star | \hat{O}_i^\dagger \hat{H} \hat{O}_j | \Psi^\star \rangle\]

where \(\hat{O}_i\) is an element of expansion_operators.

Then the generalized eigenvalue problem

\[H^{(s)} \vec{x} = E S \vec{x}\]

is solved, with \(S\) the overlap matrix:

\[S_{i,j} = \langle \Psi^\star | \hat{O}_i^\dagger \hat{O}_j | \Psi^\star \rangle\]

Parameters

• hamiltonian (SpinHamiltonian) – The Hamiltonian in its spin representation.

• state_prep_circ (Circuit) – The state preparation circuit.

• expansion_operators (list<SpinHamiltonian>) – The set of operators \({O_i}_i\) generating the subspace of interest.

• qpu (QPUHandler) – The QPU.

• nbshots (int, optional) – The number of shots. Defaults to 0: infinite number of shots.

• threshold (float, optional) – The numerical threshold.

• return_matrices (bool, optional) – If set to True, the function returns the matrices \(H^{(s)}\) and \(S\). Defaults to False.

Returns

Improved energy provided by the QSE procedure. matrix_h (Optional[np.ndarray]): The subspace Hamiltonian \(H^{(s)}\). Only if return_matrices is True. matrix_s (Optional[np.ndarray]): The overlap matrix \(S\). Only if return_matrices is True.

Return type

e_qse (float)

Example:

import numpy as np
from qat.core import Term
from qat.fermion import SpinHamiltonian
from qat.lang.AQASM import Program, RY, CNOT, RZ
from qat.qpus import get_default_qpu
from qat.plugins import SeqOptim

# We instantiate the Hamiltonian we want to approximate the ground state energy of
hamiltonian = SpinHamiltonian(2, [Term(1, op, [0, 1]) for op in ["XX", "YY", "ZZ"]])

# We construct the variational circuit (ansatz)
prog = Program()
reg = prog.qalloc(2)
theta = [prog.new_var(float, '\\theta_%s'%i) for i in range(3)]
RY(theta[0])(reg[0])
RY(theta[1])(reg[1])
RZ(theta[2])(reg[1])
CNOT(reg[0], reg[1])
circ = prog.to_circ()

# Construct a (variational) job with the variational circuit and the observable
job = circ.to_job(observable=hamiltonian,
                nbshots=0)

# We now build a stack that can handle variational jobs
qpu = get_default_qpu()

optimizer = SeqOptim(ncycles=10, x0=[0, 0.5, 0])
stack = optimizer | qpu

# We submit the job and print the optimized variational energy (the exact GS energy is -3)
result = stack.submit(job)
E_min = -3
print("E(VQE) = %s (err = %s %%)"%(result.value, 100*abs((result.value-E_min)/E_min)))
e_vqe = result.value

# We use the optimal parameters found by VQE
opt_circ = circ.bind_variables(eval(result.meta_data["parameter_map"]))

expansion_operators = [SpinHamiltonian(2, [], 1.0),
                    SpinHamiltonian(2, [Term(1., "ZZ", [0, 1])])]

from qat.fermion.chemistry.qse import apply_quantum_subspace_expansion
e_qse = apply_quantum_subspace_expansion(hamiltonian,
                                        opt_circ,
                                        expansion_operators,
                                        qpu,
                                        return_matrices=False)

print("E(QSE) = %s (err = %s %%)"%(e_qse, 100*abs((e_qse-E_min)/E_min)))

E(VQE) = -3.0 (err = 0.0 %)
E(QSE) = -3.0 (err = 0.0 %)

Unitary Coupled-Cluster (UCC)

qat.fermion.chemistry.ucc.transform_integrals_to_new_basis(one_body_integrals: ndarray, two_body_integrals: ndarray, u_mat: ndarray) → Tuple[ndarray, ndarray]

Change one and two body integrals (indices p, q…) to new basis (indices i, j…) using transformation U such that

\[\hat{c}_{i}=\sum_{q}U_{qi}c_{q}\]

i.e

\[\hat{I}_{ij} =\sum_{pq}U_{pi}I_{pq}U_{jq}^{\dagger} \hat{I}_{ijkl}=\sum_{pqrs}U_{pi}U_{qj}I_{pqrs}U_{kr}^{\dagger}U_{ls}^{\dagger}\]

Parameters

• one_body_integrals (np.ndarray) – One-body integrals \(I_{pq}\).

• two_body_integrals (np.ndarray) – Two-body integrals \(I_{pqrs}\).

• u_mat (np.ndarray) – Transformation matrix \(U\).

Returns

• h_hat_ij (np.ndarray): One-body integrals \(\hat{I}_{ij}\).

• h_hat_ijkl (np.ndarray): Two-body integrals \(\hat{I}_{ijkl}\).

Return type

Tuple[np.ndarray, np.ndarray]

qat.fermion.chemistry.ucc.compute_active_space_integrals(one_body_integrals: ndarray, two_body_integrals: ndarray, active_indices: List[int], occupied_indices: List[int]) → Tuple[ndarray, ndarray, float]

Restrict one- and two-body integrals for given list of active indices.

\[ \begin{align}\begin{aligned}\forall u,v\in \mathcal{A},\; I^{(a)}_{uv} = I_{uv} + \sum_{i\in \mathcal{O}} 2 I_{i,u,v,i} - I_{i,u,i,v}\\\forall u,v,w,x \in \mathcal{A}, I^{(a)}_{uvwx} = I_{uvwx}\\c^{(a)} = c + \sum_{i\in\mathcal{O}} I_{ii} + \sum_{ij\in\mathcal{O}} 2I_{ijji} - I_{ijij}\end{aligned}\end{align} \]

Parameters

• one_body_integrals (np.ndarray) – Array of one-body integrals \(I_{uv}\). Must be 2D.

• two_body_integrals (np.ndarray) – Array of two-body integrals \(I_{uvwx}\). Must be 4D.

• active_indices (List[int]) – Active indices.

• occupied_indices (List[int]) – Occupied indices.

Returns

• 2D array of one-body integrals \(I_{uv}^{(a)}\),

• 4D array of two-body integrals \(I_{uvwx}^{(a)}\),

• core constant \(c^{(a)}\).

Return type

Tuple[np.ndarray, np.ndarray, float]

qat.fermion.chemistry.ucc.convert_to_h_integrals(one_body_integrals: ndarray, two_body_integrals: ndarray) → Tuple[ndarray, ndarray]

Convert from \(I_{uv},I_{uvwx}\) to \(h_{pq},h_{pqrs}\), with

\[ \begin{align}\begin{aligned}h_{u\sigma, v\sigma'} = I_{u, v} \delta_{\sigma, \sigma'}\\h_{u\sigma_1, v\sigma_2, w\sigma_2', x\sigma_1'} = I_{uvwx} \left((1-\delta_{\sigma,\sigma'}) + \delta_{\sigma,\sigma'} (1-\delta_{u,v})(1-\delta_{w,x}) \right)\end{aligned}\end{align} \]

and where the one- and two-body integrals are defined as:

\[I_{uv}\equiv(u|h|v)=\int\mathrm{d}r\phi_{u}^{*}(r)T\phi_{v}(r)\]

\[I_{uvwx}\equiv(ux|vw)=\iint\mathrm{d}r_{1}\mathrm{d}r_{2}\phi_{u}^{*}(r_{1})\phi_{x}(r_{1})v(r_{12})\phi_{v}^{*}(r_{2})\phi_{w}(r_{2})\]

with \(T\) (resp. \(v\)) the one- (resp. two-) body potentials, and \(\phi_u(r)\) is the molecular orbital wavefunction.

The \(h\) integrals are used to construct hamiltonians of the ElectronicStructureHamiltonian type.

Parameters

• one_body_integrals (np.ndarray) – Array of one-body integrals \(I_{uv}\). Must be 2D.

• two_body_integrals (np.ndarray) – Array of two-body integrals \(I_{uvwx}\). Must be 4D.

Returns

• the \(h_{pq}\) integrals,

• the \(h_{pqrs}\) integrals.

Return type

Tuple[np.ndarray, np.ndarray]

qat.fermion.chemistry.ucc.construct_ucc_ansatz(cluster_ops: List[SpinHamiltonian], ket_hf: int, n_steps: int = 1) → Program

Builds the parametric state preparation circuit implementing the provided cluster operator.

The returned function maps \(\vec{\theta}\) to a QRoutine describing \(Q\) such as:

\[Q \vert \vec{0} \rangle = \vert \mathrm{UCC} (\vec{\theta}) \rangle = e^{T(\vec{\theta})} \vert \mathrm{HF}\rangle\]

Parameters

• cluster_ops (List[SpinHamiltonian]) – The cluster operators iT (note the i factor).

• ket_hf (int) – The Hartree-Fock state in integer representation.

• n_steps (int) – Number of trotter steps.

Returns

The parametric program implementing the UCCSD method.

Return type

Program

qat.fermion.chemistry.ucc.select_active_orbitals(noons: List[float], n_electrons: int, threshold_1: Optional[float] = 0.02, threshold_2: Optional[float] = 0.001) → Tuple[List[int], List[int]]

Selects the right active space and freezes core electrons according to their NOONs.

This function is an implementation of the Complete Active Space (CAS) approach. It divides orbital space into sets of active and inactive orbitals, the occupation number of the latter remaining unchanged during the computation.

Parameters

• noons (np.ndarray) – The natural orbital occupation numbers in descending order (from high occupations to low occupations)

• n_electrons (int) – The number of electrons.

• threshold_1 (Optional[float]) – The upper threshold \(\varepsilon_1\) on the NOON of an active orbital. Defaults to 0.02.

• threshold_2 (Optional[float]) – The lower threshold \(\varepsilon_2\) on the NOON of an active orbital. Defaults to 0.001.

Returns

• active_so (List[int]): The list of active spatial orbitals.

• inactive_occupied_so (List[int]): The list of core spatial orbitals.

Return type

Tuple[List[int], List[int]]

qat.fermion.chemistry.ucc.guess_init_params(two_body_integrals: ndarray, n_electrons: int, orbital_energies: List[float], noons: Optional[List[float]] = None) → List[float]

Find initial parameters using Møller-Plesset perturbation theory.

The trial parametrization is efficiently improved upon the Hartree-Fock solution (which would set every initial parameter to zero) thanks to the following formula identifying the UCC parameters in the Møller-Plesset (MP2) solution :

\[\theta_a^i = 0\]

\[\theta_{a, b}^{i, j} = \frac{h_{a, b, i, j} - h_{a, b, j, i}}{\epsilon_i + \epsilon_j -\epsilon_a - \epsilon_b}\]

where \(h_{p, q, r, s}\) is the 2-electron molecular orbital integral, and \(\epsilon_i\) is the orbital energy.

Parameters

• two_body_integrals (np.ndarray) – 4D array of two-body integrals \(I_{uvwx}\).

• n_electrons (int) – The number of active electrons of the system.

• noons (List[float]) – the natural-orbital occupation numbers \(n_i\), sorted in descending order (from high occupations to low occupations) (doubled due to spin degeneracy).

• orbital_energies (List[float]) – The energies of the molecular orbitals \(\epsilon_i\) (doubled due to spin degeneracy).

Returns

The list of initial coefficients \(\{\theta_{a}^{i}, a \in \mathcal{I}', i \in \mathcal{O}' \} \cup \{\theta_{ab}^{ij}, a>b, i>j, a,b \in \mathcal{I}', i,j \in \mathcal{O}'\}\),

Return type

theta_list (List[float])

qat.fermion.chemistry.ucc.get_hf_ket(n_electrons: int, nqbits: int) → int

Get Hartree-Fock state stored as a vector with right-to-left orbitals indexing.

Parameters

• n_electrons (int) – The number of active electrons of the system.

• nqbits – The number of qubits.

Returns

Hartree-Fock state.

Return type

int

qat.fermion.chemistry.ucc.get_cluster_ops(n_electrons: int, nqbits: Optional[int] = None, noons: Optional[List[float]] = None) → List[FermionHamiltonian]

Compute the cluster operators.

Parameters

• n_electrons (int) – The number of active electrons of the system.

• nqbits (Optional[int]) – The number of qubits.

• noons (Optional[List[float]]) – The natural-orbital occupation numbers \(n_i\), sorted in descending order (from high occupations to low occupations) (doubled due to spin degeneracy).

Returns

The list of cluster operators \(\{T_{a}^{i}, a \in \mathcal{I}', i \in \mathcal{O}' \} \cup \{T_{ab}^{ij}, a>b, i>j, a,b \in \mathcal{I}', i,j \in \mathcal{O}'\}\)

Return type

List[FermionHamiltonian]

Note

This function accepts as input the number of qubits or the noons. One of them is needed for the computation of the cluster operators. n_electrons and n_qbits must be pair.

Utility functions

qat.fermion.trotterisation.make_trotterisation_routine(hamiltonian: Union[SpinHamiltonian, FermionHamiltonian, ElectronicStructureHamiltonian], n_trotter_steps: int, final_time: Optional[float] = 1.0, method: Optional[str] = 'jordan-wigner') → QRoutine

This function first trotterizes the evolution operator \(e^{-i H t}\) of a Hamiltonian \(H\) using a first order approximation. If the Hamiltonian is fermionic, it is converted to its spin representation.

Parameters

• hamiltonian (Union[SpinHamiltonian, FermionHamiltonian, ElectronicStructureHamiltonian]) – Hamiltonian to trotterize.

• n_trotter_steps (int) – Number \(n\) of Trotter steps.

• final_time (Optional[float]) – Time \(t\) in the evolution operator.

• method (Optional[str]) – Method to use for the transformation to a spin representation. Other available methods include "bravyi-kitaev" and "parity". Defaults to "jordan-wigner".

Returns

Gates to apply to perform the time evolution of the chemical Hamiltonian with trotterisation.

Return type

QRoutine

Notes

• In the fermionic case :

  \[e^{-i H t} \approx \prod_{k=1}^{n} \left( \prod_{pq} e^{-i \frac{t}{n} h_{pq} c_p^\dagger c_q} \prod_{pqrs} e^{-\frac{i}{2}\frac{t}{n} h_{pqrs} e^{-i c_p^\dagger c_q^\dagger c_r c_s} } \right)\]

  This operator is then mapped to a product of Pauli operators via a Jordan-Wigner transformation and the resulting QRoutine is returned.

• The QRoutine implements a first order Trotter approximation, but higher order approximations are possible.

qat.fermion.chemistry.pyscf_tools.perform_pyscf_computation(geometry: list, basis: str, spin: int, charge: int, run_fci: bool = False)

Perform various calculations using PySCF. This function is a helper function meant to kickstart molecule studies. Its use is completely optional, and using other methods or packages is entirely possible.

This function will compute:

• The reduced density matrix,

• The orbital energies,

• The nuclear repulsion constant,

• The number of electrons,

• The one- and two-body integrals,

• The groundstate energies obtained through Hartree-Fock and 2nd order Möller-Plesset perturbation approach,

• (Optional) The groundstate energy using the full configuration interaction (full CI) approach.

Note

• The FCI computation is very expensive for big molecules. Enable it only for small molecules !

Parameters

• geometry (list) –

  Defines the molecular structure. The internal format is PySCF format:

  atom = [[atom1, (x, y, z)],
          [atom2, (x, y, z)],
          ...
          [atomN, (x, y, z)]]
  

• basis (str) – Defines the basis set.

• spin (int) – 2S, number of alpha electrons - number beta electrons to control multiplicity. If spin is None, multiplicity

• molecule. (will be guessed based on the neutral) –

• charge (int) – Charge of molecule. Affects the electron numbers.

• run_fci (bool, optional) – Whether the groundstates energies should also be computed using a full CI approach. Defaults to False.

Returns

• rdm1 (np.ndarray): Reduced density matrix.

• orbital_energies (list): List of orbital energies.

• nuclear_repulsion (float): Nuclear repulsion constant.

• nels (int): Number of electrons.

• one_body_integrals (np.ndarray): One-body integral.

• two_body_integrals (np.ndarray): Two-body integral.

• info (dict): Dictionary containing the Hartree-Fock and 2nd order Möller-Plesset computed ground state energies (and optionally the Full CI energy if run_fci is set to True).

Return type

Tuple[np.ndarray, list, float, int, np.ndarray, np.ndarray, dict]

Plugins

• AdaptVQEPlugin

• GradientDescentOptimizer

• SeqOptim

• MultipleLaunchesAnalyzer

• ZeroNoiseExtrapolator