Introduction to FQE

{ }

View on QuantumAI

Run in Google Colab

View source on GitHub

Download notebook

OpenFermion-FQE is an emulator for quantum computing specialized for simulations of Fermionic many-body problems, where FQE stands for 'Fermionic quantum emulator.' By focusing on Fermionic physics, OpenFermion-FQE realizes far more efficient emulation of quantum computing than generic quantum emulators such as Cirq, both in computation and memory costs; the speed-up and improved memory footprint originate from the use of the spin and number symmetries as well as highly optimized special algorithms.

The examples of the problems that can be simulated by OpenFermion-FQE include those in molecular electronic structure, condensed matter physics, and nuclear physics.

The initial version of OpenFermion-FQE has been developed in collaboration between QSimulate and Google Quantum AI. The source code is found in the GitHub repository (https://github.com/quantumlib/OpenFermion-FQE).

This tutorial will describe the data structures and conventions of the library.

try:
    import fqe
except ImportError:
    !pip install fqe --quiet

import fqe
import numpy as np

The FQE Wavefunction

The Wavefunction is an interface to the objects that hold the actual wavefunction data. As mentioned, the wavefunction is partitioned into sectors with fixed particle and \(Sz\) quantum number. This partitioning information is the necessary information for initializing a Wavefunction object.

As an example, we consider initializing a wavefunction with four spatial orbitals, four electrons, and different \(Sz\) expectation values. The Wavefunction object takes a list of these sectors [[n_electrons, sz, n_orbits]].

wfn = fqe.Wavefunction([[4, 4, 4], [4, 2, 4], [4, 0, 4], [4, -2, 4], [4, -4, 4]])

This command initializes a wavefunction with the following block structure:

Image of wf sectors

Each sector corresponds to a set of bit strings

\[ \vert I \rangle = \vert I_{\alpha}I_{\beta}\rangle \]

that encode a fixed particle number and fixed \(Sz\) expectation. The coefficients associated with the bitstrings in these sectors are formed into matrices. This helps with efficient vectorized computations. The row-index of the array corresponds to the \(\alpha\) spin-orbital number occupation index and the column-index corresponds to the \(\beta\)-strings. The Wavefunction object provides tools to access sectors or perform basic mathematical operations on this vector.

Methods to initialize wavefunctions

FQE wavefunctions can be initialized by calling the constructor directly.

wfn_fqe = fqe.Wavefunction([[2, -2, 4]], broken=None)

When wavefunctions are first created, they are initialized to empty values. We can see this by printing out the wavefunction.

wfn_fqe.print_wfn()

Sector N = 2 : S_z = -2

To set the values of a wavefunction, we can use the set_wfn method with a provided strategy.

wfn_fqe.set_wfn(strategy="hartree-fock")
wfn_fqe.print_wfn()

Sector N = 2 : S_z = -2
a'0000'b'0011' (1+0j)

Users can access the wavefunction through the get_sector method. This returns the entire matrix of data representing the specified sector of the wavefunction. For example, we can grab the sector corresponding to \(n = 2\) and \(sz = -2\) by doing the following.

interesting_states = wfn_fqe.get_coeff((2, -2))
print(interesting_states)

[[1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]]

Other than the Wavefunction constructor, several utility methods are available to initialize wavefunctions. The function fqe.get_wavefunction builds a wavefunction with definite particle number and spin.

wfn_fqe = fqe.get_wavefunction(4, -2, 10)

The function fqe.get_wavefunction_multiple constructs multiple wavefunctions with different particle number, spin, and orbital number.

wfn_fqe1, wfn_fqe2 = fqe.get_wavefunction_multiple([[4, 0, 10], [5, -5, 10]])

There are also functions like fqe.get_number_conserving_wavefunction and fqe.get_spin_conserving_wavefunction to get number or spin conserving wavefunctions, respectively.

# Get a spin conserving wavefunction.
spin_conserving_wfn = fqe.get_spin_conserving_wavefunction(2, 4)

# Get a number conserving wavefunction.
number_conserving_wfn = fqe.get_number_conserving_wavefunction(2, 4)

Conversions between FQE and Cirq wavefunction representations

Wavefunctions on \(n\) qubits in Cirq are represented by Numpy arrays with \(2^n\) amplitudes.

nqubits = 4

cirq_wfn = np.random.rand(2**nqubits) + 1.0j * np.random.rand(2**nqubits)
cirq_wfn /= np.linalg.norm(cirq_wfn)

print("Cirq wavefunction:")
print(*cirq_wfn, sep="\n")

Cirq wavefunction:
(0.013896449806861794+0.23221431680562818j)
(0.166469937176125+0.1677724510037618j)
(0.09075944422730575+0.23578905933003944j)
(0.22166559274324935+0.2757340874479102j)
(0.06664512169623714+0.28076381705891745j)
(0.006604516899126683+0.2264991185104466j)
(0.1944380035556502+0.18799406075328381j)
(0.09531075263861329+0.09913398939488889j)
(0.19435205690407295+0.07208206745571455j)
(0.13177959280983956+0.030091830199165205j)
(0.2882523823401827+0.13235227097707394j)
(0.15710789387178647+0.2965705083877339j)
(0.10972180327683587+0.0549056883676476j)
(0.13463044259659268+0.03915330817748353j)
(0.04268126178644539+0.2620894369245064j)
(0.26746752372622873+0.15199737633910376j)

To convert from this representation to the FQE representation, the function fqe.from_cirq can be used.

fqe_wfn = fqe.from_cirq(cirq_wfn, thresh=0.0001)
fqe_wfn.print_wfn()

Sector N = 0 : S_z = 0
a'00'b'00' (0.013896449806861794+0.23221431680562818j)
Sector N = 1 : S_z = -1
a'00'b'01' (0.06664512169623714+0.28076381705891745j)
a'00'b'10' (0.166469937176125+0.1677724510037618j)
Sector N = 1 : S_z = 1
a'01'b'00' (0.19435205690407295+0.07208206745571455j)
a'10'b'00' (0.09075944422730575+0.23578905933003944j)
Sector N = 2 : S_z = -2
a'00'b'11' (0.006604516899126683+0.2264991185104466j)
Sector N = 2 : S_z = 0
a'01'b'01' (0.10972180327683587+0.0549056883676476j)
a'01'b'10' (0.13177959280983956+0.030091830199165205j)
a'10'b'01' (-0.1944380035556502-0.18799406075328381j)
a'10'b'10' (0.22166559274324935+0.2757340874479102j)
Sector N = 2 : S_z = 2
a'11'b'00' (0.2882523823401827+0.13235227097707394j)
Sector N = 3 : S_z = -1
a'01'b'11' (0.13463044259659268+0.03915330817748353j)
a'10'b'11' (-0.09531075263861329-0.09913398939488889j)
Sector N = 3 : S_z = 1
a'11'b'01' (-0.04268126178644539-0.2620894369245064j)
a'11'b'10' (0.15710789387178647+0.2965705083877339j)
Sector N = 4 : S_z = 0
a'11'b'11' (-0.26746752372622873-0.15199737633910376j)

Note: The thresh argument is the value below which amplitudes are considered zero.

We can convert back to the Cirq representation using fqe._to_cirq.

cirq_wfn_from_fqe = fqe.to_cirq(fqe_wfn)

print("Cirq wavefunction from FQE:")
print(*cirq_wfn_from_fqe, sep="\n")

Cirq wavefunction from FQE:
(0.013896449806861794+0.23221431680562818j)
(0.166469937176125+0.1677724510037618j)
(0.09075944422730575+0.23578905933003944j)
(0.22166559274324935+0.2757340874479102j)
(0.06664512169623714+0.28076381705891745j)
(0.006604516899126683+0.2264991185104466j)
(0.1944380035556502+0.18799406075328381j)
(0.09531075263861329+0.09913398939488889j)
(0.19435205690407295+0.07208206745571455j)
(0.13177959280983956+0.030091830199165205j)
(0.2882523823401827+0.13235227097707394j)
(0.15710789387178647+0.2965705083877339j)
(0.10972180327683587+0.0549056883676476j)
(0.13463044259659268+0.03915330817748353j)
(0.04268126178644539+0.2620894369245064j)
(0.26746752372622873+0.15199737633910376j)

assert np.allclose(cirq_wfn_from_fqe, cirq_wfn)

An important thing to note in these conversions is the ordering of the \(\alpha\) and \(\beta\) strings in the converted wavefunctions. The FQE uses the OpenFermion convention of interleaved \(\alpha\) and \(\beta\) orbitals. Thus when converting to Cirq we first convert each bitstring into an OpenFermion operator and then call normal ordering.

Printing and saving wavefunctions

Printing is currently available as alpha beta strings followed by the coefficient as well as orbital occupation representation.

print('String formatting')
fqe_wfn.print_wfn(fmt='str')

String formatting
Sector N = 0 : S_z = 0
a'00'b'00' (0.013896449806861794+0.23221431680562818j)
Sector N = 1 : S_z = -1
a'00'b'01' (0.06664512169623714+0.28076381705891745j)
a'00'b'10' (0.166469937176125+0.1677724510037618j)
Sector N = 1 : S_z = 1
a'01'b'00' (0.19435205690407295+0.07208206745571455j)
a'10'b'00' (0.09075944422730575+0.23578905933003944j)
Sector N = 2 : S_z = -2
a'00'b'11' (0.006604516899126683+0.2264991185104466j)
Sector N = 2 : S_z = 0
a'01'b'01' (0.10972180327683587+0.0549056883676476j)
a'01'b'10' (0.13177959280983956+0.030091830199165205j)
a'10'b'01' (-0.1944380035556502-0.18799406075328381j)
a'10'b'10' (0.22166559274324935+0.2757340874479102j)
Sector N = 2 : S_z = 2
a'11'b'00' (0.2882523823401827+0.13235227097707394j)
Sector N = 3 : S_z = -1
a'01'b'11' (0.13463044259659268+0.03915330817748353j)
a'10'b'11' (-0.09531075263861329-0.09913398939488889j)
Sector N = 3 : S_z = 1
a'11'b'01' (-0.04268126178644539-0.2620894369245064j)
a'11'b'10' (0.15710789387178647+0.2965705083877339j)
Sector N = 4 : S_z = 0
a'11'b'11' (-0.26746752372622873-0.15199737633910376j)

print('Occupation formatting')
fqe_wfn.print_wfn(fmt='occ')

Occupation formatting
Sector N = 0 : S_z = 0
.. (0.013896449806861794+0.23221431680562818j)
Sector N = 1 : S_z = -1
.b (0.06664512169623714+0.28076381705891745j)
b. (0.166469937176125+0.1677724510037618j)
Sector N = 1 : S_z = 1
.a (0.19435205690407295+0.07208206745571455j)
a. (0.09075944422730575+0.23578905933003944j)
Sector N = 2 : S_z = -2
bb (0.006604516899126683+0.2264991185104466j)
Sector N = 2 : S_z = 0
.2 (0.10972180327683587+0.0549056883676476j)
ba (0.13177959280983956+0.030091830199165205j)
ab (-0.1944380035556502-0.18799406075328381j)

2. (0.22166559274324935+0.2757340874479102j)
Sector N = 2 : S_z = 2
aa (0.2882523823401827+0.13235227097707394j)
Sector N = 3 : S_z = -1
b2 (0.13463044259659268+0.03915330817748353j)
2b (-0.09531075263861329-0.09913398939488889j)
Sector N = 3 : S_z = 1
a2 (-0.04268126178644539-0.2620894369245064j)
2a (0.15710789387178647+0.2965705083877339j)
Sector N = 4 : S_z = 0
22 (-0.26746752372622873-0.15199737633910376j)

Wavefunctions can also be saved to disk using the save method which takes a filename and optional path.

Action on Wavefunctions: Fermionic algebra operations and their unitaries on the state

FermionOperators can be directly passed in to create a new wavefunction based on application of the operators. The FermionOperators are parsed according to the interleaved \(\alpha\) \(\beta\) indexing of the spin-orbitals. This means that odd index FermionOperators correspond to \(\beta\)-spin orbitals and even are \(\alpha\)-spin orbitals.

Sharp Edge:

The user must be careful to not break the symmetry of the wavefunction. If a request to apply an operator to a state takes the wavefunction outside of the specified symmetry sector the FQE will not execute the command. Effectively, the FQE requires the user to have more knowledge of what type of operations their Wavefunction object can support.

from openfermion import FermionOperator, hermitian_conjugated

ops = FermionOperator('2^ 0', 1.2)

new_wfn = fqe_wfn.apply(ops + hermitian_conjugated(ops))

Unitary operations

Any simulator backend must be able to perform unitary evolution on a state. The FQE accomplishes this by implementing code for evolving a state via the action of a unitary generated by fermionic generators. Given a fermion operator \(g\), the unitary

\[ e^{-i (g + g^{\dagger})} \]

can be applied to the state. It can be shown that this evolution operator can be rewritten as

\[ e^{-i(g + g^{\dagger})\epsilon } = \cos\left(\epsilon\right) \mathbb{I}_{s}(gg^{\dagger}) + \cos\left(\epsilon\right) \mathbb{I}_{s}(g^{\dagger}g) - i\sin\left(\epsilon\right) \left(g + g^{\dagger}\right) \left[\mathbb{I}_{s}(gg^{\dagger}) + \mathbb{I}_{s}(g^{\dagger}g)\right] + \mathbb{I}_{!s} \]

The \(\mathbb{I}_{!s}\) is for setting the coefficients of the unitary that are not in the subspace \(\mathcal{H}_{s} \subset \mathcal{H}\) where \(gg^{\dagger}\) is 0.

The user can specify a fermionic monomial in OpenFermion and use the time_evolve method of the Wavefunction object to call the evolution. All the rules for preserving symmetries must be maintained as before.

i, j, theta = 0, 1, np.pi / 3
op = (FermionOperator(((2 * i, 1), (2 * j, 0)), coefficient=-1j * theta) + 
      FermionOperator(((2 * j, 1), (2 * i, 0)), coefficient=1j * theta))

new_wfn = fqe_wfn.time_evolve(1.0, op)

new_wfn.print_wfn()

Sector N = 0 : S_z = 0
a'00'b'00' (0.013896449806861794+0.23221431680562818j)
Sector N = 1 : S_z = -1
a'00'b'01' (0.06664512169623714+0.28076381705891745j)
a'00'b'10' (0.166469937176125+0.1677724510037618j)
Sector N = 1 : S_z = 1
a'01'b'00' (0.0185760441178328-0.16815828158639307j)
a'10'b'00' (0.21369354067033883+0.1803194312389721j)
Sector N = 2 : S_z = -2
a'00'b'11' (0.006604516899126683+0.2264991185104466j)
Sector N = 2 : S_z = 0
a'01'b'01' (0.22324915217874+0.1902604765567627j)
a'01'b'10' (-0.12607823805566964-0.22374680931962754j)
a'10'b'01' (-0.0021971327910466076-0.04644730943798737j)
a'10'b'10' (0.2249572714453149+0.16392733312279995j)
Sector N = 2 : S_z = 2
a'11'b'00' (0.2882523823401827+0.13235227097707394j)
Sector N = 3 : S_z = -1
a'01'b'11' (0.14985675433715018+0.10542920728321267j)
a'10'b'11' (0.0689380070920852-0.015659235173542722j)
Sector N = 3 : S_z = 1
a'11'b'01' (-0.04268126178644539-0.2620894369245064j)
a'11'b'10' (0.15710789387178647+0.2965705083877339j)
Sector N = 4 : S_z = 0
a'11'b'11' (-0.26746752372622873-0.15199737633910376j)

In other tutorials we will do a deeper dive into supported time-evolution operations. To serve a full functioning time-evolution platform the FQE also implements arbitrary time-evolution of full Hamiltonian operators consisting of sums of non-commuting terms. The Taylor and Chebyshev expansion methods are used to do the exact time evolution.