 View on QuantumAI
|
 Run in Google Colab
|
 View source on GitHub
|
 Download notebook
|
Setup
Install the OpenFermion package:
try:
import openfermion
except ImportError:
!pip install git+https://github.com/quantumlib/OpenFermion.git@master#egg=openfermion
Initializing the FermionOperator data structure
Fermionic systems are often treated in second quantization where arbitrary operators can be expressed using the fermionic creation and annihilation operators, \(a^\dagger_k\) and \(a_k\). The fermionic ladder operators play a similar role to their qubit ladder operator counterparts, \(\sigma^+_k\) and \(\sigma^-_k\) but are distinguished by the canonical fermionic anticommutation relations, \(\{a^\dagger_i, a^\dagger_j\} = \{a_i, a_j\} = 0\) and \(\{a_i, a_j^\dagger\} = \delta_{ij}\). Any weighted sums of products of these operators are represented with the FermionOperator data structure in OpenFermion. The following are examples of valid FermionOperators:
\[ \begin{align} & a_1 \nonumber \\ & 1.7 a^\dagger_3 \nonumber \\ &-1.7 \, a^\dagger_3 a_1 \nonumber \\ &(1 + 2i) \, a^\dagger_4 a^\dagger_3 a_9 a_1 \nonumber \\ &(1 + 2i) \, a^\dagger_4 a^\dagger_3 a_9 a_1 - 1.7 \, a^\dagger_3 a_1 \nonumber \end{align} \]
The FermionOperator class is contained in \(\textrm{ops/_fermion_operator.py}\). In order to support fast addition of FermionOperator instances, the class is implemented as hash table (python dictionary). The keys of the dictionary encode the strings of ladder operators and values of the dictionary store the coefficients. The strings of ladder operators are encoded as a tuple of 2-tuples which we refer to as the "terms tuple". Each ladder operator is represented by a 2-tuple. The first element of the 2-tuple is an int indicating the tensor factor on which the ladder operator acts. The second element of the 2-tuple is Boole: 1 represents raising and 0 represents lowering. For instance, \(a^\dagger_8\) is represented in a 2-tuple as \((8, 1)\). Note that indices start at 0 and the identity operator is an empty list. Below we give some examples of operators and their terms tuple:
\[ \begin{align} I & \mapsto () \nonumber \\ a_1 & \mapsto ((1, 0),) \nonumber \\ a^\dagger_3 & \mapsto ((3, 1),) \nonumber \\ a^\dagger_3 a_1 & \mapsto ((3, 1), (1, 0)) \nonumber \\ a^\dagger_4 a^\dagger_3 a_9 a_1 & \mapsto ((4, 1), (3, 1), (9, 0), (1, 0)) \nonumber \end{align} \]
Note that when initializing a single ladder operator one should be careful to add the comma after the inner pair. This is because in python ((1, 2)) = (1, 2) whereas ((1, 2),) = ((1, 2),). The "terms tuple" is usually convenient when one wishes to initialize a term as part of a coded routine. However, the terms tuple is not particularly intuitive. Accordingly, OpenFermion also supports another user-friendly, string notation below. This representation is rendered when calling "print" on a FermionOperator.
\[ \begin{align} I & \mapsto \textrm{""} \nonumber \\ a_1 & \mapsto \textrm{"1"} \nonumber \\ a^\dagger_3 & \mapsto \textrm{"3^"} \nonumber \\ a^\dagger_3 a_1 & \mapsto \textrm{"3^}\;\textrm{1"} \nonumber \\ a^\dagger_4 a^\dagger_3 a_9 a_1 & \mapsto \textrm{"4^}\;\textrm{3^}\;\textrm{9}\;\textrm{1"} \nonumber \end{align} \]
Let's initialize our first term! We do it two different ways below.
from openfermion.ops import FermionOperator
my_term = FermionOperator(((3, 1), (1, 0)))
print(my_term)
my_term = FermionOperator('3^ 1')
print(my_term)
1.0 [3^ 1] 1.0 [3^ 1]
The preferred way to specify the coefficient in openfermion is to provide an optional coefficient argument. If not provided, the coefficient defaults to 1. In the code below, the first method is preferred. The multiplication in the second method actually creates a copy of the term, which introduces some additional cost. All inplace operands (such as +=) modify classes whereas binary operands such as + create copies. Important caveats are that the empty tuple FermionOperator(()) and the empty string FermionOperator('') initializes identity. The empty initializer FermionOperator() initializes the zero operator.
good_way_to_initialize = FermionOperator('3^ 1', -1.7)
print(good_way_to_initialize)
bad_way_to_initialize = -1.7 * FermionOperator('3^ 1')
print(bad_way_to_initialize)
identity = FermionOperator('')
print(identity)
zero_operator = FermionOperator()
print(zero_operator)
-1.7 [3^ 1] -1.7 [3^ 1] 1.0 [] 0
Note that FermionOperator has only one attribute: .terms. This attribute is the dictionary which stores the term tuples.
my_operator = FermionOperator('4^ 1^ 3 9', 1. + 2.j)
print(my_operator)
print(my_operator.terms)
(1+2j) [4^ 1^ 3 9]
{((4, 1), (1, 1), (3, 0), (9, 0)): (1+2j)}
Manipulating the FermionOperator data structure
So far we have explained how to initialize a single FermionOperator such as \(-1.7 \, a^\dagger_3 a_1\). However, in general we will want to represent sums of these operators such as \((1 + 2i) \, a^\dagger_4 a^\dagger_3 a_9 a_1 - 1.7 \, a^\dagger_3 a_1\). To do this, just add together two FermionOperators! We demonstrate below.
from openfermion.ops import FermionOperator
term_1 = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator = term_1 + term_2
print(my_operator)
my_operator = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator += term_2
print('')
print(my_operator)
-1.7 [3^ 1] + (1+2j) [4^ 3^ 9 1] -1.7 [3^ 1] + (1+2j) [4^ 3^ 9 1]
The print function prints each term in the operator on a different line. Note that the line my_operator = term_1 + term_2 creates a new object, which involves a copy of term_1 and term_2. The second block of code uses the inplace method +=, which is more efficient. This is especially important when trying to construct a very large FermionOperator. FermionOperators also support a wide range of builtins including, str(), repr(), ==, !=, *=, *, /, /=, +, +=, -, -=, - and **. Note that since FermionOperators involve floats, == and != check for (in)equality up to numerical precision. We demonstrate some of these methods below.
term_1 = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator = term_1 - 33. * term_2
print(my_operator)
my_operator *= 3.17 * (term_2 + term_1) ** 2
print('')
print(my_operator)
print('')
print(term_2 ** 3)
print('')
print(term_1 == 2.*term_1 - term_1)
print(term_1 == my_operator)
56.1 [3^ 1] + (1+2j) [4^ 3^ 9 1] 513.9489299999999 [3^ 1 3^ 1 3^ 1] + (-302.32289999999995-604.6457999999999j) [3^ 1 3^ 1 4^ 3^ 9 1] + (-302.32289999999995-604.6457999999999j) [3^ 1 4^ 3^ 9 1 3^ 1] + (-533.511+711.348j) [3^ 1 4^ 3^ 9 1 4^ 3^ 9 1] + (9.161299999999999+18.322599999999998j) [4^ 3^ 9 1 3^ 1 3^ 1] + (16.166999999999998-21.555999999999997j) [4^ 3^ 9 1 3^ 1 4^ 3^ 9 1] + (16.166999999999998-21.555999999999997j) [4^ 3^ 9 1 4^ 3^ 9 1 3^ 1] + (-34.87-6.34j) [4^ 3^ 9 1 4^ 3^ 9 1 4^ 3^ 9 1] -4.912999999999999 [3^ 1 3^ 1 3^ 1] True False
Additionally, there are a variety of methods that act on the FermionOperator data structure. We demonstrate a small subset of those methods here.
from openfermion.utils import commutator, count_qubits, hermitian_conjugated
from openfermion.transforms import normal_ordered
# Get the Hermitian conjugate of a FermionOperator, count its qubit, check if it is normal-ordered.
term_1 = FermionOperator('4^ 3 3^', 1. + 2.j)
print(hermitian_conjugated(term_1))
print(term_1.is_normal_ordered())
print(count_qubits(term_1))
# Normal order the term.
term_2 = normal_ordered(term_1)
print('')
print(term_2)
print(term_2.is_normal_ordered())
# Compute a commutator of the terms.
print('')
print(commutator(term_1, term_2))
(1-2j) [3 3^ 4] False 5 (1+2j) [4^] + (-1-2j) [4^ 3^ 3] True (-3+4j) [4^ 3 3^ 4^] + (3-4j) [4^ 3 3^ 4^ 3^ 3] + (-3+4j) [4^ 3^ 3 4^ 3 3^] + (3-4j) [4^ 4^ 3 3^]
The QubitOperator data structure
The QubitOperator data structure is another essential part of openfermion. As the name suggests, QubitOperator is used to store qubit operators in almost exactly the same way that FermionOperator is used to store fermion operators. For instance \(X_0 Z_3 Y_4\) is a QubitOperator. The internal representation of this as a terms tuple would be \(((0, \textrm{"X"}), (3, \textrm{"Z"}), (4, \textrm{"Y"}))\). Note that one important difference between QubitOperator and FermionOperator is that the terms in QubitOperator are always sorted in order of tensor factor. In some cases, this enables faster manipulation. We initialize some QubitOperators below.
from openfermion.ops import QubitOperator
my_first_qubit_operator = QubitOperator('X1 Y2 Z3')
print(my_first_qubit_operator)
print(my_first_qubit_operator.terms)
operator_2 = QubitOperator('X3 Z4', 3.17)
operator_2 -= 77. * my_first_qubit_operator
print('')
print(operator_2)
1.0 [X1 Y2 Z3]
{((1, 'X'), (2, 'Y'), (3, 'Z')): 1.0}
-77.0 [X1 Y2 Z3] +
3.17 [X3 Z4]
Jordan-Wigner and Bravyi-Kitaev
openfermion provides functions for mapping FermionOperators to QubitOperators.
from openfermion.ops import FermionOperator
from openfermion.transforms import jordan_wigner, bravyi_kitaev
from openfermion.utils import hermitian_conjugated
from openfermion.linalg import eigenspectrum
# Initialize an operator.
fermion_operator = FermionOperator('2^ 0', 3.17)
fermion_operator += hermitian_conjugated(fermion_operator)
print(fermion_operator)
# Transform to qubits under the Jordan-Wigner transformation and print its spectrum.
jw_operator = jordan_wigner(fermion_operator)
print('')
print(jw_operator)
jw_spectrum = eigenspectrum(jw_operator)
print(jw_spectrum)
# Transform to qubits under the Bravyi-Kitaev transformation and print its spectrum.
bk_operator = bravyi_kitaev(fermion_operator)
print('')
print(bk_operator)
bk_spectrum = eigenspectrum(bk_operator)
print(bk_spectrum)
3.17 [0^ 2] + 3.17 [2^ 0] (1.585+0j) [X0 Z1 X2] + (1.585+0j) [Y0 Z1 Y2] [-3.17 -3.17 0. 0. 0. 0. 3.17 3.17] (1.585+0j) [X0 Y1 Y2] + (-1.585+0j) [Y0 Y1 X2] [-3.17 -3.17 0. 0. 0. 0. 3.17 3.17]
We see that despite the different representation, these operators are iso-spectral. We can also apply the Jordan-Wigner transform in reverse to map arbitrary QubitOperators to FermionOperators. Note that we also demonstrate the .compress() method (a method on both FermionOperators and QubitOperators) which removes zero entries.
from openfermion.transforms import reverse_jordan_wigner
# Initialize QubitOperator.
my_operator = QubitOperator('X0 Y1 Z2', 88.)
my_operator += QubitOperator('Z1 Z4', 3.17)
print(my_operator)
# Map QubitOperator to a FermionOperator.
mapped_operator = reverse_jordan_wigner(my_operator)
print('')
print(mapped_operator)
# Map the operator back to qubits and make sure it is the same.
back_to_normal = jordan_wigner(mapped_operator)
back_to_normal.compress()
print('')
print(back_to_normal)
88.0 [X0 Y1 Z2] + 3.17 [Z1 Z4] 3.17 [] + -88j [1 0] + 88j [1 0^] + 88j [1^ 0] + -88j [1^ 0^] + -6.34 [1^ 1] + 176j [2^ 2 1 0] + -176j [2^ 2 1 0^] + -176j [2^ 2 1^ 0] + 176j [2^ 2 1^ 0^] + -6.34 [4^ 4] + 12.68 [4^ 4 1^ 1] 88.0 [X0 Y1 Z2] + 3.17 [Z1 Z4]
Sparse matrices and the Hubbard model
Often, one would like to obtain a sparse matrix representation of an operator which can be analyzed numerically. There is code in both openfermion.transforms and openfermion.utils which facilitates this. The function get_sparse_operator converts either a FermionOperator, a QubitOperator or other more advanced classes such as InteractionOperator to a scipy.sparse.csc matrix. There are numerous functions in openfermion.utils which one can call on the sparse operators such as "get_gap", "get_hartree_fock_state", "get_ground_state", etc. We show this off by computing the ground state energy of the Hubbard model. To do that, we use code from the openfermion.hamiltonians module which constructs lattice models of fermions such as Hubbard models.
from openfermion.hamiltonians import fermi_hubbard
from openfermion.linalg import get_sparse_operator, get_ground_state
from openfermion.transforms import jordan_wigner
# Set model.
x_dimension = 2
y_dimension = 2
tunneling = 2.
coulomb = 1.
magnetic_field = 0.5
chemical_potential = 0.25
periodic = 1
spinless = 1
# Get fermion operator.
hubbard_model = fermi_hubbard(
x_dimension, y_dimension, tunneling, coulomb, chemical_potential,
magnetic_field, periodic, spinless)
print(hubbard_model)
# Get qubit operator under Jordan-Wigner.
jw_hamiltonian = jordan_wigner(hubbard_model)
jw_hamiltonian.compress()
print('')
print(jw_hamiltonian)
# Get scipy.sparse.csc representation.
sparse_operator = get_sparse_operator(hubbard_model)
print('')
print(sparse_operator)
print('\nEnergy of the model is {} in units of T and J.'.format(
get_ground_state(sparse_operator)[0]))
-0.25 [0^ 0] +
1.0 [0^ 0 1^ 1] +
1.0 [0^ 0 2^ 2] +
-2.0 [0^ 1] +
-2.0 [0^ 2] +
-2.0 [1^ 0] +
-0.25 [1^ 1] +
1.0 [1^ 1 3^ 3] +
-2.0 [1^ 3] +
-2.0 [2^ 0] +
-0.25 [2^ 2] +
1.0 [2^ 2 3^ 3] +
-2.0 [2^ 3] +
-2.0 [3^ 1] +
-2.0 [3^ 2] +
-0.25 [3^ 3]
0.5 [] +
-1.0 [X0 X1] +
-1.0 [X0 Z1 X2] +
-1.0 [Y0 Y1] +
-1.0 [Y0 Z1 Y2] +
-0.375 [Z0] +
0.25 [Z0 Z1] +
0.25 [Z0 Z2] +
-1.0 [X1 Z2 X3] +
-1.0 [Y1 Z2 Y3] +
-0.375 [Z1] +
0.25 [Z1 Z3] +
-1.0 [X2 X3] +
-1.0 [Y2 Y3] +
-0.375 [Z2] +
0.25 [Z2 Z3] +
-0.375 [Z3]
<Compressed Sparse Column sparse matrix of dtype 'complex128'
with 47 stored elements and shape (16, 16)>
Coords Values
(1, 1) (-0.25+0j)
(2, 1) (-2+0j)
(4, 1) (-2+0j)
(1, 2) (-2+0j)
(2, 2) (-0.25+0j)
(8, 2) (-2+0j)
(3, 3) (0.5+0j)
(6, 3) (2+0j)
(9, 3) (-2+0j)
(1, 4) (-2+0j)
(4, 4) (-0.25+0j)
(8, 4) (-2+0j)
(5, 5) (0.5+0j)
(6, 5) (-2+0j)
(9, 5) (-2+0j)
(3, 6) (2+0j)
(5, 6) (-2+0j)
(6, 6) (-0.5+0j)
(10, 6) (-2+0j)
(12, 6) (2+0j)
(7, 7) (1.25+0j)
(11, 7) (-2+0j)
(13, 7) (2+0j)
(2, 8) (-2+0j)
(4, 8) (-2+0j)
(8, 8) (-0.25+0j)
(3, 9) (-2+0j)
(5, 9) (-2+0j)
(9, 9) (-0.5+0j)
(10, 9) (-2+0j)
(12, 9) (-2+0j)
(6, 10) (-2+0j)
(9, 10) (-2+0j)
(10, 10) (0.5+0j)
(7, 11) (-2+0j)
(11, 11) (1.25+0j)
(14, 11) (2+0j)
(6, 12) (2+0j)
(9, 12) (-2+0j)
(12, 12) (0.5+0j)
(7, 13) (2+0j)
(13, 13) (1.25+0j)
(14, 13) (-2+0j)
(11, 14) (2+0j)
(13, 14) (-2+0j)
(14, 14) (1.25+0j)
(15, 15) (3+0j)
Energy of the model is -4.249999999999997 in units of T and J.
Hamiltonians in the plane wave basis
A user can write plugins to openfermion which allow for the use of, e.g., third-party electronic structure package to compute molecular orbitals, Hamiltonians, energies, reduced density matrices, coupled cluster amplitudes, etc using Gaussian basis sets. We may provide scripts which interface between such packages and openfermion in future but do not discuss them in this tutorial.
When using simpler basis sets such as plane waves, these packages are not needed. openfermion comes with code which computes Hamiltonians in the plane wave basis. Note that when using plane waves, one is working with the periodized Coulomb operator, best suited for condensed phase calculations such as studying the electronic structure of a solid. To obtain these Hamiltonians one must choose to study the system without a spin degree of freedom (spinless), one must the specify dimension in which the calculation is performed (n_dimensions, usually 3), one must specify how many plane waves are in each dimension (grid_length) and one must specify the length scale of the plane wave harmonics in each dimension (length_scale) and also the locations and charges of the nuclei. One can generate these models with plane_wave_hamiltonian() found in openfermion.hamiltonians. For simplicity, below we compute the Hamiltonian in the case of zero external charge (corresponding to the uniform electron gas, aka jellium). We also demonstrate that one can transform the plane wave Hamiltonian using a Fourier transform without effecting the spectrum of the operator.
from openfermion.hamiltonians import jellium_model
from openfermion.utils import Grid
from openfermion.linalg import eigenspectrum
from openfermion.transforms import jordan_wigner, fourier_transform
# Let's look at a very small model of jellium in 1D.
grid = Grid(dimensions=1, length=3, scale=1.0)
spinless = True
# Get the momentum Hamiltonian.
momentum_hamiltonian = jellium_model(grid, spinless)
momentum_qubit_operator = jordan_wigner(momentum_hamiltonian)
momentum_qubit_operator.compress()
print(momentum_qubit_operator)
# Fourier transform the Hamiltonian to the position basis.
position_hamiltonian = fourier_transform(momentum_hamiltonian, grid, spinless)
position_qubit_operator = jordan_wigner(position_hamiltonian)
position_qubit_operator.compress()
print('')
print (position_qubit_operator)
# Check the spectra to make sure these representations are iso-spectral.
spectral_difference = eigenspectrum(momentum_qubit_operator) - eigenspectrum(position_qubit_operator)
print('')
print(spectral_difference)
19.50047638754088 [] + -9.71044945799746 [Z0] + -0.07957747154594767 [Z0 Z1] + -0.07957747154594767 [Z0 Z2] + 0.15915494309189535 [Z1] + -0.07957747154594767 [Z1 Z2] + -9.71044945799746 [Z2] 19.500476387540854 [] + -3.289868133696451 [X0 X1] + -3.289868133696454 [X0 Z1 X2] + -3.289868133696451 [Y0 Y1] + -3.289868133696454 [Y0 Z1 Y2] + -6.420581324301009 [Z0] + -0.07957747154594766 [Z0 Z1] + -0.07957747154594763 [Z0 Z2] + -3.289868133696451 [X1 X2] + -3.289868133696451 [Y1 Y2] + -6.4205813243010095 [Z1] + -0.07957747154594766 [Z1 Z2] + -6.420581324301009 [Z2] [2.75474088e-14 2.66175970e-14 2.84217094e-14 7.10542736e-15 2.84217094e-14 2.13162821e-14 1.42108547e-14 7.10542736e-15]
Basics of MolecularData class
Data from electronic structure calculations can be saved in an OpenFermion data structure called MolecularData, which makes it easy to access within our library. Often, one would like to analyze a chemical series or look at many different Hamiltonians and sometimes the electronic structure calculations are either expensive to compute or difficult to converge (e.g. one needs to mess around with different types of SCF routines to make things converge). Accordingly, we anticipate that users will want some way to automatically database the results of their electronic structure calculations so that important data (such as the SCF integrals) can be looked up on-the-fly if the user has computed them in the past. OpenFermion supports a data provenance strategy which saves key results of the electronic structure calculation (including pointers to files containing large amounts of data, such as the molecular integrals) in an HDF5 container.
The MolecularData class stores information about molecules. One initializes a MolecularData object by specifying parameters of a molecule such as its geometry, basis, multiplicity, charge and an optional string describing it. One can also initialize MolecularData simply by providing a string giving a filename where a previous MolecularData object was saved in an HDF5 container. One can save a MolecularData instance by calling the class's .save() method. This automatically saves the instance in a data folder specified during OpenFermion installation. The name of the file is generated automatically from the instance attributes and optionally provided description. Alternatively, a filename can also be provided as an optional input if one wishes to manually name the file.
When electronic structure calculations are run, the data files for the molecule can be automatically updated. If one wishes to later use that data they either initialize MolecularData with the instance filename or initialize the instance and then later call the .load() method.
Basis functions are provided to initialization using a string such as "6-31g". Geometries can be specified using a simple txt input file (see geometry_from_file function in molecular_data.py) or can be passed using a simple python list format demonstrated below. Atoms are specified using a string for their atomic symbol. Distances should be provided in angstrom. Below we initialize a simple instance of MolecularData without performing any electronic structure calculations.
from openfermion.chem import MolecularData
# Set parameters to make a simple molecule.
diatomic_bond_length = .7414
geometry = [('H', (0., 0., 0.)), ('H', (0., 0., diatomic_bond_length))]
basis = 'sto-3g'
multiplicity = 1
charge = 0
description = str(diatomic_bond_length)
# Make molecule and print out a few interesting facts about it.
molecule = MolecularData(geometry, basis, multiplicity,
charge, description)
print('Molecule has automatically generated name {}'.format(
molecule.name))
print('Information about this molecule would be saved at:\n{}\n'.format(
molecule.filename))
print('This molecule has {} atoms and {} electrons.'.format(
molecule.n_atoms, molecule.n_electrons))
for atom, atomic_number in zip(molecule.atoms, molecule.protons):
print('Contains {} atom, which has {} protons.'.format(
atom, atomic_number))
Molecule has automatically generated name H2_sto-3g_singlet_0.7414 Information about this molecule would be saved at: /tmpfs/src/tf_docs_env/lib/python3.10/site-packages/openfermion/testing/data/H2_sto-3g_singlet_0.7414 This molecule has 2 atoms and 2 electrons. Contains H atom, which has 1 protons. Contains H atom, which has 1 protons.
If we had previously computed this molecule using an electronic structure package, we can call molecule.load() to populate all sorts of interesting fields in the data structure. Though we make no assumptions about what electronic structure packages users might install, we assume that the calculations are saved in OpenFermion's MolecularData objects. Currently plugins are available for Psi4 (OpenFermion-Psi4) and PySCF (OpenFermion-PySCF), and there may be more in the future. For the purposes of this example, we will load data that ships with OpenFermion to make a plot of the energy surface of hydrogen. Note that helper functions to initialize some interesting chemical benchmarks are found in openfermion.utils.
# Set molecule parameters.
basis = 'sto-3g'
multiplicity = 1
bond_length_interval = 0.1
n_points = 25
# Generate molecule at different bond lengths.
hf_energies = []
fci_energies = []
bond_lengths = []
for point in range(3, n_points + 1):
bond_length = bond_length_interval * point
bond_lengths += [bond_length]
description = str(round(bond_length,2))
print(description)
geometry = [('H', (0., 0., 0.)), ('H', (0., 0., bond_length))]
molecule = MolecularData(
geometry, basis, multiplicity, description=description)
# Load data.
molecule.load()
# Print out some results of calculation.
print('\nAt bond length of {} angstrom, molecular hydrogen has:'.format(
bond_length))
print('Hartree-Fock energy of {} Hartree.'.format(molecule.hf_energy))
print('MP2 energy of {} Hartree.'.format(molecule.mp2_energy))
print('FCI energy of {} Hartree.'.format(molecule.fci_energy))
print('Nuclear repulsion energy between protons is {} Hartree.'.format(
molecule.nuclear_repulsion))
for orbital in range(molecule.n_orbitals):
print('Spatial orbital {} has energy of {} Hartree.'.format(
orbital, molecule.orbital_energies[orbital]))
hf_energies += [molecule.hf_energy]
fci_energies += [molecule.fci_energy]
# Plot.
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(0)
plt.plot(bond_lengths, fci_energies, 'x-')
plt.plot(bond_lengths, hf_energies, 'o-')
plt.ylabel('Energy in Hartree')
plt.xlabel('Bond length in angstrom')
plt.show()
0.3 At bond length of 0.30000000000000004 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.59382776458494 Hartree. MP2 energy of -0.5997818888874376 Hartree. FCI energy of -0.6018037168352988 Hartree. Nuclear repulsion energy between protons is 1.7639240286333335 Hartree. Spatial orbital 0 has energy of -0.8028666187118976 Hartree. Spatial orbital 1 has energy of 1.368938795250223 Hartree. 0.4 At bond length of 0.4 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.904361397713845 Hartree. MP2 energy of -0.9114367296997896 Hartree. FCI energy of -0.9141497082141279 Hartree. Nuclear repulsion energy between protons is 1.322943021475 Hartree. Spatial orbital 0 has energy of -0.7452125332909346 Hartree. Spatial orbital 1 has energy of 1.167416395038123 Hartree. 0.5 At bond length of 0.5 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.0429962765072203 Hartree. MP2 energy of -1.0514858618835026 Hartree. FCI energy of -1.055159796496619 Hartree. Nuclear repulsion energy between protons is 1.05835441718 Hartree. Spatial orbital 0 has energy of -0.6908223275124182 Hartree. Spatial orbital 1 has energy of 0.9886736693741557 Hartree. 0.6 At bond length of 0.6000000000000001 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.1011282431368283 Hartree. MP2 energy of -1.1113317674794059 Hartree. FCI energy of -1.1162860078265333 Hartree. Nuclear repulsion energy between protons is 0.8819620143166668 Hartree. Spatial orbital 0 has energy of -0.6408762643991083 Hartree. Spatial orbital 1 has energy of 0.8380849781488554 Hartree. 0.7 At bond length of 0.7000000000000001 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.1173490350703141 Hartree. MP2 energy of -1.1295815518324286 Hartree. FCI energy of -1.1361894542708848 Hartree. Nuclear repulsion energy between protons is 0.7559674408428572 Hartree. Spatial orbital 0 has energy of -0.5954634716745077 Hartree. Spatial orbital 1 has energy of 0.7141652810054977 Hartree. 0.8 At bond length of 0.8 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.1108503969854677 Hartree. MP2 energy of -1.125453091983669 Hartree. FCI energy of -1.1341476663578605 Hartree. Nuclear repulsion energy between protons is 0.6614715107375 Hartree. Spatial orbital 0 has energy of -0.5544958797639515 Hartree. Spatial orbital 1 has energy of 0.6126180834930599 Hartree. 0.9 At bond length of 0.9 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.0919140401071208 Hartree. MP2 energy of -1.1092713614102743 Hartree. FCI energy of -1.1205602806186359 Hartree. Nuclear repulsion energy between protons is 0.5879746762111111 Hartree. Spatial orbital 0 has energy of -0.5176680309627193 Hartree. Spatial orbital 1 has energy of 0.5284772431612719 Hartree. 1.0 At bond length of 1.0 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.0661086480817525 Hartree. MP2 energy of -1.0866648035682114 Hartree. FCI energy of -1.101150329303588 Hartree. Nuclear repulsion energy between protons is 0.52917720859 Hartree. Spatial orbital 0 has energy of -0.4844416789615343 Hartree. Spatial orbital 1 has energy of 0.4575019361641323 Hartree. 1.1 At bond length of 1.1 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.0365388735395633 Hartree. MP2 energy of -1.0608057828548 Hartree. FCI energy of -1.0791929438802303 Hartree. Nuclear repulsion energy between protons is 0.4810701896272727 Hartree. Spatial orbital 0 has energy of -0.4542186927189922 Hartree. Spatial orbital 1 has energy of 0.39669590841169317 Hartree. 1.2 At bond length of 1.2000000000000002 angstrom, molecular hydrogen has: Hartree-Fock energy of -1.0051067048832343 Hartree. MP2 energy of -1.0336620966065275 Hartree. FCI energy of -1.056740745132591 Hartree. Nuclear repulsion energy between protons is 0.4409810071583334 Hartree. Spatial orbital 0 has energy of -0.42650264072324545 Hartree. Spatial orbital 1 has energy of 0.3441268787843519 Hartree. 1.3 At bond length of 1.3 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.9731106139491602 Hartree. MP2 energy of -1.006590197243068 Hartree. FCI energy of -1.0351862652459827 Hartree. Nuclear repulsion energy between protons is 0.40705939122307694 Hartree. Spatial orbital 0 has energy of -0.40094651614872323 Hartree. Spatial orbital 1 has energy of 0.29851388983893856 Hartree. 1.4 At bond length of 1.4000000000000001 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.9414806527840186 Hartree. MP2 energy of -0.9805687328589369 Hartree. FCI energy of -1.0154682481426915 Hartree. Nuclear repulsion energy between protons is 0.3779837204214286 Hartree. Spatial orbital 0 has energy of -0.3773228239693589 Hartree. Spatial orbital 1 has energy of 0.2589019723128326 Hartree. 1.5 At bond length of 1.5 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.9108735526174462 Hartree. MP2 energy of -0.9562942016363287 Hartree. FCI energy of -0.9981493524136993 Hartree. Nuclear repulsion energy between protons is 0.35278480572666665 Hartree. Spatial orbital 0 has energy of -0.3554774880192603 Hartree. Spatial orbital 1 has energy of 0.22449543728518556 Hartree. 1.6 At bond length of 1.6 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.8817324479521221 Hartree. MP2 energy of -0.9342411695719502 Hartree. FCI energy of -0.9834727280932716 Hartree. Nuclear repulsion energy between protons is 0.33073575536875 Hartree. Spatial orbital 0 has energy of -0.3352963505891839 Hartree. Spatial orbital 1 has energy of 0.19459795539978458 Hartree. 1.7 At bond length of 1.7000000000000002 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.8543376249706873 Hartree. MP2 energy of -0.9147126228432917 Hartree. FCI energy of -0.9714266876511379 Hartree. Nuclear repulsion energy between protons is 0.31128071093529414 Hartree. Spatial orbital 0 has energy of -0.3166860441012991 Hartree. Spatial orbital 1 has energy of 0.16859732072886116 Hartree. 1.8 At bond length of 1.8 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.8288481459846262 Hartree. MP2 energy of -0.8978805380405488 Hartree. FCI energy of -0.9618169521199167 Hartree. Nuclear repulsion energy between protons is 0.29398733810555555 Hartree. Spatial orbital 0 has energy of -0.29956322022346127 Hartree. Spatial orbital 1 has energy of 0.1459602966246219 Hartree. 1.9 At bond length of 1.9000000000000001 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.8053328430090216 Hartree. MP2 energy of -0.8838158086555916 Hartree. FCI energy of -0.9543388534525794 Hartree. Nuclear repulsion energy between protons is 0.2785143203105263 Hartree. Spatial orbital 0 has energy of -0.2838478652114843 Hartree. Spatial orbital 1 has energy of 0.1262260489187021 Hartree. 2.0 At bond length of 2.0 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.783792652466311 Hartree. MP2 energy of -0.8725100604810707 Hartree. FCI energy of -0.9486411117424077 Hartree. Nuclear repulsion energy between protons is 0.264588604295 Hartree. Spatial orbital 0 has energy of -0.26945922244444087 Hartree. Spatial orbital 1 has energy of 0.1089973681318946 Hartree. 2.1 At bond length of 2.1 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.7641776498896591 Hartree. MP2 energy of -0.8638929333599986 Hartree. FCI energy of -0.944374680781336 Hartree. Nuclear repulsion energy between protons is 0.25198914694761904 Hartree. Spatial orbital 0 has energy of -0.25631405863438766 Hartree. Spatial orbital 1 has energy of 0.0939315946347687 Hartree. 2.2 At bond length of 2.2 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.7464013483548488 Hartree. MP2 energy of -0.8578472435835268 Hartree. FCI energy of -0.9412240334330747 Hartree. Nuclear repulsion energy between protons is 0.24053509481363636 Hartree. Spatial orbital 0 has energy of -0.24432699142489445 Hartree. Spatial orbital 1 has energy of 0.08073270789777327 Hartree. 2.3 At bond length of 2.3000000000000003 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.7303533198131085 Hartree. MP2 energy of -0.854222721280043 Hartree. FCI energy of -0.938922385789505 Hartree. Nuclear repulsion energy between protons is 0.2300770472130435 Hartree. Spatial orbital 0 has energy of -0.2334122089960839 Hartree. Spatial orbital 1 has energy of 0.069144997403965 Hartree. 2.4 At bond length of 2.4000000000000004 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.7159100590078187 Hartree. MP2 energy of -0.8528478930177863 Hartree. FCI energy of -0.9372549528606704 Hartree. Nuclear repulsion energy between protons is 0.2204905035791667 Hartree. Spatial orbital 0 has energy of -0.22348569714218677 Hartree. Spatial orbital 1 has energy of 0.05894803199690316 Hartree. 2.5 At bond length of 2.5 angstrom, molecular hydrogen has: Hartree-Fock energy of -0.7029435983728299 Hartree. MP2 energy of -0.8535394695489013 Hartree. FCI energy of -0.9360549198442759 Hartree. Nuclear repulsion energy between protons is 0.21167088343600005 Hartree. Spatial orbital 0 has energy of -0.2144671917818725 Hartree. Spatial orbital 1 has energy of 0.04995242289757026 Hartree.
The geometry data needed to generate MolecularData can also be retrieved from the PubChem online database by inputting the molecule's name.
from openfermion.chem import geometry_from_pubchem
methane_geometry = geometry_from_pubchem('methane')
print(methane_geometry)
[('C', (0, 0, 0)), ('H', (0.5541, 0.7996, 0.4965)), ('H', (0.6833, -0.8134, -0.2536)), ('H', (-0.7782, -0.3735, 0.6692)), ('H', (-0.4593, 0.3874, -0.9121))]
InteractionOperator and InteractionRDM for efficient numerical representations
Fermion Hamiltonians can be expressed as \(H = h_0 + \sum_{pq} h_{pq}\, a^\dagger_p a_q + \frac{1}{2} \sum_{pqrs} h_{pqrs} \, a^\dagger_p a^\dagger_q a_r a_s\) where \(h_0\) is a constant shift due to the nuclear repulsion and \(h_{pq}\) and \(h_{pqrs}\) are the famous molecular integrals. Since fermions interact pairwise, their energy is thus a unique function of the one-particle and two-particle reduced density matrices which are expressed in second quantization as \(\rho_{pq} = \left \langle p \mid a^\dagger_p a_q \mid q \right \rangle\) and \(\rho_{pqrs} = \left \langle pq \mid a^\dagger_p a^\dagger_q a_r a_s \mid rs \right \rangle\), respectively.
Because the RDMs and molecular Hamiltonians are both compactly represented and manipulated as 2- and 4- index tensors, we can represent them in a particularly efficient form using similar data structures. The InteractionOperator data structure can be initialized for a Hamiltonian by passing the constant \(h_0\) (or 0), as well as numpy arrays representing \(h_{pq}\) (or \(\rho_{pq}\)) and \(h_{pqrs}\) (or \(\rho_{pqrs}\)). Importantly, InteractionOperators can also be obtained by calling MolecularData.get_molecular_hamiltonian() or by calling the function get_interaction_operator() (found in openfermion.transforms) on a FermionOperator. The InteractionRDM data structure is similar but represents RDMs. For instance, one can get a molecular RDM by calling MolecularData.get_molecular_rdm(). When generating Hamiltonians from the MolecularData class, one can choose to restrict the system to an active space.
These classes inherit from the same base class, PolynomialTensor. This data structure overloads the slice operator [] so that one can get or set the key attributes of the InteractionOperator: \(\textrm{.constant}\), \(\textrm{.one_body_coefficients}\) and \(\textrm{.two_body_coefficients}\) . For instance, InteractionOperator[(p, 1), (q, 1), (r, 0), (s, 0)] would return \(h_{pqrs}\) and InteractionRDM would return \(\rho_{pqrs}\). Importantly, the class supports fast basis transformations using the method PolynomialTensor.rotate_basis(rotation_matrix). But perhaps most importantly, one can map the InteractionOperator to any of the other data structures we've described here.
Below, we load MolecularData from a saved calculation of LiH. We then obtain an InteractionOperator representation of this system in an active space. We then map that operator to qubits. We then demonstrate that one can rotate the orbital basis of the InteractionOperator using random angles to obtain a totally different operator that is still iso-spectral.
from openfermion.chem import MolecularData
from openfermion.transforms import get_fermion_operator, jordan_wigner
from openfermion.linalg import get_ground_state, get_sparse_operator
import numpy
import scipy
import scipy.linalg
# Load saved file for LiH.
diatomic_bond_length = 1.45
geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., diatomic_bond_length))]
basis = 'sto-3g'
multiplicity = 1
# Set Hamiltonian parameters.
active_space_start = 1
active_space_stop = 3
# Generate and populate instance of MolecularData.
molecule = MolecularData(geometry, basis, multiplicity, description="1.45")
molecule.load()
# Get the Hamiltonian in an active space.
molecular_hamiltonian = molecule.get_molecular_hamiltonian(
occupied_indices=range(active_space_start),
active_indices=range(active_space_start, active_space_stop))
# Map operator to fermions and qubits.
fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian)
qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)
qubit_hamiltonian.compress()
print('The Jordan-Wigner Hamiltonian in canonical basis follows:\n{}'.format(qubit_hamiltonian))
# Get sparse operator and ground state energy.
sparse_hamiltonian = get_sparse_operator(qubit_hamiltonian)
energy, state = get_ground_state(sparse_hamiltonian)
print('Ground state energy before rotation is {} Hartree.\n'.format(energy))
# Randomly rotate.
n_orbitals = molecular_hamiltonian.n_qubits // 2
n_variables = int(n_orbitals * (n_orbitals - 1) / 2)
numpy.random.seed(1)
random_angles = numpy.pi * (1. - 2. * numpy.random.rand(n_variables))
kappa = numpy.zeros((n_orbitals, n_orbitals))
index = 0
for p in range(n_orbitals):
for q in range(p + 1, n_orbitals):
kappa[p, q] = random_angles[index]
kappa[q, p] = -numpy.conjugate(random_angles[index])
index += 1
# Build the unitary rotation matrix.
difference_matrix = kappa + kappa.transpose()
rotation_matrix = scipy.linalg.expm(kappa)
# Apply the unitary.
molecular_hamiltonian.rotate_basis(rotation_matrix)
# Get qubit Hamiltonian in rotated basis.
qubit_hamiltonian = jordan_wigner(molecular_hamiltonian)
qubit_hamiltonian.compress()
print('The Jordan-Wigner Hamiltonian in rotated basis follows:\n{}'.format(qubit_hamiltonian))
# Get sparse Hamiltonian and energy in rotated basis.
sparse_hamiltonian = get_sparse_operator(qubit_hamiltonian)
energy, state = get_ground_state(sparse_hamiltonian)
print('Ground state energy after rotation is {} Hartree.'.format(energy))
The Jordan-Wigner Hamiltonian in canonical basis follows: -7.49894690201071 [] + -0.0029329964409502266 [X0 X1 Y2 Y3] + 0.0029329964409502266 [X0 Y1 Y2 X3] + 0.01291078027311749 [X0 Z1 X2] + -0.0013743761078958677 [X0 Z1 X2 Z3] + 0.011536413200774975 [X0 X2] + 0.0029329964409502266 [Y0 X1 X2 Y3] + -0.0029329964409502266 [Y0 Y1 X2 X3] + 0.01291078027311749 [Y0 Z1 Y2] + -0.0013743761078958677 [Y0 Z1 Y2 Z3] + 0.011536413200774975 [Y0 Y2] + 0.16199475388004184 [Z0] + 0.011536413200774975 [Z0 X1 Z2 X3] + 0.011536413200774975 [Z0 Y1 Z2 Y3] + 0.12444770133137588 [Z0 Z1] + 0.054130445793298836 [Z0 Z2] + 0.05706344223424907 [Z0 Z3] + 0.012910780273117487 [X1 Z2 X3] + -0.0013743761078958677 [X1 X3] + 0.012910780273117487 [Y1 Z2 Y3] + -0.0013743761078958677 [Y1 Y3] + 0.16199475388004186 [Z1] + 0.05706344223424907 [Z1 Z2] + 0.054130445793298836 [Z1 Z3] + -0.013243698330265966 [Z2] + 0.08479609543670981 [Z2 Z3] + -0.013243698330265952 [Z3] Ground state energy before rotation is -7.862773163027987 Hartree. The Jordan-Wigner Hamiltonian in rotated basis follows: -7.498946902010706 [] + -0.024260054468246924 [X0 X1 Y2 Y3] + 0.024260054468246924 [X0 Y1 Y2 X3] + -0.08262397170394281 [X0 Z1 X2] + -0.01673498917401365 [X0 Z1 X2 Z3] + -0.005524724636333846 [X0 X2] + 0.024260054468246924 [Y0 X1 X2 Y3] + -0.024260054468246924 [Y0 Y1 X2 X3] + -0.08262397170394281 [Y0 Z1 Y2] + -0.01673498917401365 [Y0 Z1 Y2 Z3] + -0.005524724636333846 [Y0 Y2] + 0.042483580038933266 [Z0] + -0.005524724636333846 [Z0 X1 Z2 X3] + -0.005524724636333846 [Z0 Y1 Z2 Y3] + 0.08238127510732211 [Z0 Z1] + 0.05413044579329884 [Z0 Z2] + 0.07839050026154576 [Z0 Z3] + -0.08262397170394281 [X1 Z2 X3] + -0.01673498917401365 [X1 X3] + -0.08262397170394281 [Y1 Z2 Y3] + -0.01673498917401365 [Y1 Y3] + 0.042483580038933266 [Z1] + 0.07839050026154576 [Z1 Z2] + 0.05413044579329884 [Z1 Z3] + 0.1062674755108426 [Z2] + 0.0842084056061702 [Z2 Z3] + 0.1062674755108426 [Z3] Ground state energy after rotation is -7.862773163027982 Hartree.
Quadratic Hamiltonians and Slater determinants
The general electronic structure Hamiltonian
\(H = h_0 + \sum_{pq} h_{pq}\, a^\dagger_p a_q + \frac{1}{2} \sum_{pqrs} h_{pqrs} \, a^\dagger_p a^\dagger_q a_r a_s\) contains terms that act on up to 4 sites, or
is quartic in the fermionic creation and annihilation operators. However, in many situations
we may fruitfully approximate these Hamiltonians by replacing these quartic terms with
terms that act on at most 2 fermionic sites, or quadratic terms, as in mean-field approximation theory.
These Hamiltonians have a number of
special properties one can exploit for efficient simulation and manipulation of the Hamiltonian, thus
warranting a special data structure. We refer to Hamiltonians which
only contain terms that are quadratic in the fermionic creation and annihilation operators
as quadratic Hamiltonians, and include the general case of non-particle conserving terms as in
a general Bogoliubov transformation. Eigenstates of quadratic Hamiltonians can be prepared
efficiently on both a quantum and classical computer, making them amenable to initial guesses for
many more challenging problems.
A general quadratic Hamiltonian takes the form
\[H = \sum_{p, q} (M_{pq} - \mu \delta_{pq}) a^\dagger_p a_q + \frac{1}{2} \sum_{p, q} (\Delta_{pq} a^\dagger_p a^\dagger_q + \Delta_{pq}^* a_q a_p) + \text{constant},\]
where \(M\) is a Hermitian matrix, \(\Delta\) is an antisymmetric matrix, \(\delta_{pq}\) is the Kronecker delta symbol, and \(\mu\) is a chemical potential term which we keep separate from \(M\) so that we can use it to adjust the expectation of the total number of particles. In OpenFermion, quadratic Hamiltonians are conveniently represented and manipulated using the QuadraticHamiltonian class, which stores \(M\), \(\Delta\), \(\mu\) and the constant. It is specialized to exploit the properties unique to quadratic Hamiltonians. Like InteractionOperator and InteractionRDM, it inherits from the PolynomialTensor class.
The BCS mean-field model of superconductivity is a quadratic Hamiltonian. The following code constructs an instance of this model as a FermionOperator, converts it to a QuadraticHamiltonian, and then computes its ground energy:
from openfermion.hamiltonians import mean_field_dwave
from openfermion.transforms import get_quadratic_hamiltonian
# Set model.
x_dimension = 2
y_dimension = 2
tunneling = 2.
sc_gap = 1.
periodic = True
# Get FermionOperator.
mean_field_model = mean_field_dwave(
x_dimension, y_dimension, tunneling, sc_gap, periodic=periodic)
# Convert to QuadraticHamiltonian
quadratic_hamiltonian = get_quadratic_hamiltonian(mean_field_model)
# Compute the ground energy
ground_energy = quadratic_hamiltonian.ground_energy()
print(ground_energy)
-10.000000000000005
Any quadratic Hamiltonian may be rewritten in the form
\[H = \sum_p \varepsilon_p b^\dagger_p b_p + \text{constant},\]
where the \(b_p\) are new annihilation operators that satisfy the fermionic anticommutation relations, and which are linear combinations of the old creation and annihilation operators. This form of \(H\) makes it easy to deduce its eigenvalues; they are sums of subsets of the \(\varepsilon_p\), which we call the orbital energies of \(H\). The following code computes the orbital energies and the constant:
orbital_energies, constant = quadratic_hamiltonian.orbital_energies()
print(orbital_energies)
print()
print(constant)
[1. 1. 1. 1. 4. 4. 4. 4.] -10.000000000000005
Eigenstates of quadratic hamiltonians are also known as fermionic Gaussian states, and they can be prepared efficiently on a quantum computer. One can use OpenFermion to obtain circuits for preparing these states. The following code obtains the description of a circuit which prepares the ground state (operations that can be performed in parallel are grouped together), along with a description of the starting state to which the circuit should be applied:
from openfermion.circuits import gaussian_state_preparation_circuit
circuit_description, start_orbitals = gaussian_state_preparation_circuit(quadratic_hamiltonian)
for parallel_ops in circuit_description:
print(parallel_ops)
print('')
print(start_orbitals)
('pht',)
((6, 7, 1.5707963267948966, 0.0),)
('pht', (5, 6, 1.5707963267948966, 0.0))
((4, 5, 1.0471975511965983, -3.141592653589793), (6, 7, 1.047197551196598, 3.141592653589793))
('pht', (3, 4, 1.5707963267948966, 0.0), (5, 6, 1.5707963267948966, 0.0))
((2, 3, 1.2309594173407747, -2.275957200481571e-15), (4, 5, 1.230959417340775, -7.771561172376097e-16), (6, 7, 1.1071487177940909, 3.141592653589792))
('pht', (1, 2, 1.5707963267948966, 0.0), (3, 4, 1.5707963267948966, 0.0), (5, 6, 1.5707963267948966, 0.0))
((0, 1, 1.047197551196598, -3.1415926535897927), (2, 3, 1.0471975511965979, 3.1415926535897927), (4, 5, 1.3181160716528173, -3.141592653589791), (6, 7, 1.3181160716528173, -3.608224830031759e-16))
('pht', (1, 2, 1.5707963267948966, 0.0), (3, 4, 1.5707963267948966, 0.0), (5, 6, 1.5707963267948966, 0.0))
((2, 3, 0.9553166181245089, 4.440892098500626e-16), (4, 5, 0.9553166181245089, 2.7755575615628914e-16), (6, 7, 1.10714871779409, -8.881784197001252e-16))
('pht', (3, 4, 1.5707963267948966, 0.0), (5, 6, 1.5707963267948966, 0.0))
((4, 5, 0.7853981633974482, 3.1415926535897922), (6, 7, 0.7853981633974481, 2.220446049250313e-16))
((5, 6, 1.5707963267948966, 0.0),)
[]
In the circuit description, each elementary operation is either a tuple of the form \((i, j, \theta, \varphi)\), indicating the operation \(\exp[i \varphi a_j^\dagger a_j]\exp[\theta (a_i^\dagger a_j - a_j^\dagger a_i)]\), which is a Givens rotation of modes \(i\) and \(j\), or the string 'pht', indicating the particle-hole transformation on the last fermionic mode, which is the operator \(\mathcal{B}\) such that \(\mathcal{B} a_N \mathcal{B}^\dagger = a_N^\dagger\) and leaves the rest of the ladder operators unchanged. Operations that can be performed in parallel are grouped together.
In the special case that a quadratic Hamiltonian conserves particle number (\(\Delta = 0\)), its eigenstates take the form
\[\lvert \Psi_S \rangle = b^\dagger_{1}\cdots b^\dagger_{N_f}\lvert \text{vac} \rangle,\qquad b^\dagger_{p} = \sum_{k=1}^N Q_{pq}a^\dagger_q,\]
where \(Q\) is an \(N_f \times N\) matrix with orthonormal rows. These states are also known as Slater determinants. OpenFermion also provides functionality to obtain circuits for preparing Slater determinants starting with the matrix \(Q\) as the input.