This is the first of several tutorial explaining key concepts in OpenFermion-Cirq. These tutorials build on one another and should be studied in order. In this tutorial we will discuss the compilation of circuits for implementing arbitrary rotations of the single-particle basis of an electronic structure simulation. As an example, we show how one can use these methods to simulate the evolution of an arbitrary non-interacting fermion model.
In order to represent fermionic systems on a quantum computer one must first discretize space. Usually, one expands the many-body wavefunction in a basis of spin-orbitals $\varphi_p = \varphi_p(r)$ which are single-particle basis functions. For reasons of spatial efficiency, all NISQ (and even most error-corrected) algorithms for simulating fermionic systems focus on representing operators in second-quantization. Second-quantized operators are expressed using the fermionic creation and annihilation operators, $a^\dagger_p$ and $a_p$. The action of $a^\dagger_p$ is to excite a fermion in spin-orbital $\varphi_p$ and the action of $a_p$ is to annihilate a fermion from spin-orbital $\varphi_p$. Specifically, if electron $i$ is represented in a space of spin-orbitals $\{\varphi_p(r_i)\}$ then $a^\dagger_p$ and $a_p$ are related to Slater determinants through the equivalence, $$ \langle r_0 \cdots r_{\eta-1} | a^\dagger_{0} \cdots a^\dagger_{\eta-1} | \varnothing\rangle \equiv \sqrt{\frac{1}{\eta!}} \begin{vmatrix} \varphi_{0}\left(r_0\right) & \varphi_{1}\left( r_0\right) & \cdots & \varphi_{\eta-1} \left( r_0\right) \\ \varphi_{0}\left(r_1\right) & \varphi_{1}\left( r_1\right) & \cdots & \varphi_{\eta-1} \left( r_1\right) \\ \vdots & \vdots & \ddots & \vdots\\ \varphi_{0}\left(r_{\eta-1}\right) & \varphi_{1}\left(r_{\eta-1}\right) & \cdots & \varphi_{\eta-1} \left(r_{\eta-1}\right) \end{vmatrix} $$ where $\eta$ is the number of electrons in the system, $|\varnothing \rangle$ is the Fermi vacuum and $\varphi_p(r)=\langle r|\varphi_p \rangle$ are the single-particle orbitals that define the basis. By using a basis of Slater determinants, we ensure antisymmetry in the encoded state.
Very often in electronic structure calculations one would like to rotate the single-particle basis. That is, one would like to generate new orbitals that are formed from a linear combination of the old orbitals. Any particle-conserving rotation of the single-particle basis can be expressed as $$ \tilde{\varphi}_p = \sum_{q} \varphi_q u_{pq} \quad \tilde{a}^\dagger_p = \sum_{q} a^\dagger_q u_{pq} \quad \tilde{a}_p = \sum_{q} a_q u_{pq}^* $$ where $\tilde{\varphi}_p$, $\tilde{a}^\dagger_p$, and $\tilde{a}^\dagger_p$ correspond to spin-orbitals and operators in the rotated basis and $u$ is an $N\times N$ unitary matrix. From the Thouless theorem, this single-particle rotation is equivalent to applying the $2^N \times 2^N$ operator $$ U(u) = \exp\left(\sum_{pq} \left[\log u \right]_{pq} \left(a^\dagger_p a_q - a^\dagger_q a_p\right)\right) $$ where $\left[\log u\right]_{pq}$ is the $(p, q)$ element of the matrix $\log u$.
There are many reasons that one might be interested in performing such basis rotations. For instance, one might be interested in preparing the Hartree-Fock (mean-field) state of a chemical system, by rotating from some initial orbitals (e.g. atomic orbitals or plane waves) into the molecular orbitals of the system. Alternatively, one might be interested in rotating from a basis where certain operators are diagonal (e.g. the kinetic operator is diagonal in the plane wave basis) to a basis where certain other operators are diagonal (e.g. the Coulomb operator is diagonal in the position basis). Thus, it is a very useful thing to be able to apply circuits corresponding to $U(u)$ on a quantum computer in low depth.
OpenFermion-Cirq prominently features routines for implementing the linear depth / linear connectivity basis transformations described in Phys. Rev. Lett. 120, 110501. While we will not discuss this functionality here, we also support routines for compiling the more general form of these transformations which do not conserve particle-number, known as a Bogoliubov transfomation, using routines described in Phys. Rev. Applied 9, 044036. We will not discuss the details of how these methods are implemented here and instead refer readers to those papers. All that one needs in order to compile the circuit $U(u)$ using OpenFermion-Cirq is the $N \times N$ matrix $u$, which we refer to in documentation as the "basis_transformation_matrix". Note that if one intends to apply this matrix to a computational basis state with only $\eta$ electrons, then one can reduce the number of gates required by instead supplying the $\eta \times N$ rectangular matrix that characterizes the rotation of the occupied orbitals only. OpenFermion-Cirq will automatically take advantage of this symmetry.
In this example will show how basis transforms can be used to implement exact evolution under a random Hermitian one-body fermionic operator \begin{equation} H = \sum_{pq} T_{pq} a^\dagger_p a_q. \end{equation} That is, we will compile a circuit to implement $e^{-i H t}$ for some time $t$. Of course, this is a tractable problem classically but we discuss it here since it is often useful as a subroutine for more complex quantum simulations. To accomplish this evolution, we will use basis transformations. Suppose that $u$ is the basis transformation matrix that diagonalizes $T$. Then, we could implement $e^{-i H t}$ by implementing $U(u)^\dagger (\prod_{k} e^{-i \lambda_k Z_k}) U(u)$ where $\lambda_k$ are the eigenvalues of $T$.
Below, we initialize the T matrix characterizing $H$ and then obtain the eigenvalues $\lambda_k$ and eigenvectors $u_k$ of $T$. We print out the OpenFermion FermionOperator representation of $T$.
In [1]:
import openfermion
import numpy
# Set the number of qubits in our example.
n_qubits = 3
simulation_time = 1.
random_seed = 8317
# Generate the random one-body operator.
T = openfermion.random_hermitian_matrix(n_qubits, seed=random_seed)
# Diagonalize T and obtain basis transformation matrix (aka "u").
eigenvalues, eigenvectors = numpy.linalg.eigh(T)
basis_transformation_matrix = eigenvectors.transpose()
# Print out familiar OpenFermion "FermionOperator" form of H.
H = openfermion.FermionOperator()
for p in range(n_qubits):
for q in range(n_qubits):
term = ((p, 1), (q, 0))
H += openfermion.FermionOperator(term, T[p, q])
print(H)
Now we're ready to make a circuit! First we will use OpenFermion-Cirq to generate the basis transform $U(u)$ from the basis transformation matrix $u$ by calling the Bogoliubov transform function (named as such because this function can also handle non-particle conserving basis transformations). Then, we'll apply local $Z$ rotations to phase by the eigenvalues, then we'll apply the inverse transformation. That will finish the circuit. We're just going to print out the first rotation to keep things easy-to-read, but feel free to play around with the notebook.
In [2]:
import openfermioncirq
import cirq
# Initialize the qubit register.
qubits = cirq.LineQubit.range(n_qubits)
# Start circuit with the inverse basis rotation, print out this step.
inverse_basis_rotation = cirq.inverse(openfermioncirq.bogoliubov_transform(qubits, basis_transformation_matrix))
circuit = cirq.Circuit(inverse_basis_rotation)
print(circuit)
# Add diagonal phase rotations to circuit.
for k, eigenvalue in enumerate(eigenvalues):
phase = -eigenvalue * simulation_time
circuit.append(cirq.Rz(rads=phase).on(qubits[k]))
# Finally, restore basis.
basis_rotation = openfermioncirq.bogoliubov_transform(qubits, basis_transformation_matrix)
circuit.append(basis_rotation)
Finally, we can check whether our circuit applied to a random initial state with the exact result. Print out the fidelity with the exact result.
In [3]:
# Initialize a random initial state.
initial_state = openfermion.haar_random_vector(
2 ** n_qubits, random_seed).astype(numpy.complex64)
# Numerically compute the correct circuit output.
import scipy
hamiltonian_sparse = openfermion.get_sparse_operator(H)
exact_state = scipy.sparse.linalg.expm_multiply(
-1j * simulation_time * hamiltonian_sparse, initial_state)
# Use Cirq simulator to apply circuit.
simulator = cirq.Simulator()
result = simulator.simulate(circuit, qubit_order=qubits,
initial_state=initial_state)
simulated_state = result.final_state
# Print final fidelity.
fidelity = abs(numpy.dot(simulated_state, numpy.conjugate(exact_state)))**2
print(fidelity)
Thus, we see that the circuit correctly effects the intended evolution. We can now use Cirq's compiler to output the circuit using gates native to near-term devices, and then optimize those circuits. We'll output in QASM 2.0 just to demonstrate that functionality.
In [4]:
xmon_circuit = cirq.google.optimized_for_xmon(circuit)
print(xmon_circuit.to_qasm())