Spaces, Operators, and States

The use of quantum computers for tackling problems in quantum chemistry involves a wide variety of quantum mechanical structures. Many of these are shared with quantum chemistry problems studied with classical computers. For example, the electronic Hamiltonian may be considered a second-quantized operator acting on a fermionic Hilbert space. However, quantum computing approaches often require objects and formalisms that are atypical to “conventional” quantum chemistry – for instance, fermionic operators and states must be mapped to operators acting on and states within a qubit Hilbert space.

InQuanto provides options for representing each of these objects. Broadly, three core types of objects are used - InQuanto operator and state classes represent operators and states within a given Hilbert space. Space objects (chiefly FermionSpace and QubitSpace) can be used to describe specific Hilbert spaces and then generate these states and operators. The relationship between these objects is depicted in Fig. 24.

Fig. 24 A schematic of the vector space, state and operator classes in InQuanto. Spaces generate operators and states; fermion-to-qubit mappings transform fermionic states and operators to qubit states and operators.

We first discuss fermionic and qubit spaces and how operators and states can be mapped between the two. An option for parafermionic [10] spaces is also provided by the ParaFermionSpace object. This mapping is described in the Iterative Qubit-Excitation Based VQE algorithm section. We also consider integral operators - classes for the efficient manipulation of molecular orbital integral tensors. Finally, we consider orbital transformers and optimization – several InQuanto classes for the convenient manipulation of molecular orbitals bases.

Fermionic Spaces

Given a system of \(N\) electrons described by an orthonormal basis of \(Q\) molecular spin-orbitals, or momentum space modes, where \(Q \geq N\), one can construct a fermionic Fock space [5], or “fermion space”. Each basis vector of this abstract vector space corresponds to an occupation number (ON) vector

(67)\[|\boldsymbol{\eta}\rangle = |\eta_0, \eta_1, ..., \eta_P, ..., \eta_{Q-1} \rangle\]

where an occupation number \(\eta_P = 1\) (\(0\)) if spin-orbital \(P\) is occupied (unoccupied). Each ON vector represents a canonically ordered \(N\)-electron Slater determinant. A fermion space is an abstract linear vector space equipped with the usual properties of an inner product and resolution of identity, for a given set of spin-orbitals or modes. An arbitrary vector in this space corresponds to a linear combination of ON vectors

(68)\[|\boldsymbol{\psi}\rangle = \sum_{\boldsymbol{\eta}} \psi_{\boldsymbol{\eta}} |\boldsymbol{\eta}\rangle\]

if the spin-orbitals or modes used to construct the ON vectors are orthonormal, such that the ON vectors obey the relation

(69)\[\langle\boldsymbol{\eta}'|\boldsymbol{\eta}\rangle = \prod_{P=0}^{Q-1} \delta_{\eta'_P, \eta_P}\]

This is a useful feature of fermion spaces: a well defined but zero overlap between ON vectors with different values of \(N = \sum_{P=0}^{Q-1} \eta_P\) and \(N' = \sum_{P=0}^{Q-1} \eta'_P\) allows for a consistent treatment of systems with different numbers of electrons. Thus the ON vectors form an orthonormal basis in the \(2^Q\)-dimensional fermion space, which in principle can be decomposed into a direct sum of fermionic subspaces \(\mathcal{F}(Q, N)\), where each subspace represents a different number of electrons distributed over the \(Q\) spin-orbitals.

(70)\[\mathcal{F}(Q) = \bigoplus_{N=0}^Q \mathcal{F}(Q, N)\]

The dimension of each \(\mathcal{F}(Q, N)\) subspace is equal to the number of terms in the FCI expansion \(\binom{Q}{N}\), such that the exact (FCI) wavefunction for \(N\) electrons (in a given basis) can be expressed as a vector in \(\mathcal{F}(Q, N)\).

Anticommuting creation and annihilation operators act on the fermion space as linear maps between vectors within the space. From these operators, one- and two-body interactions can be defined, along with the fermionic ON vectors, and thus the entire fermionic problem can be described once the \(N\)-electron subspace is specified. In InQuanto, this logic is followed in the FermionSpace class which represents the fermionic Fock space \(\mathcal{F}(Q)\), and contains related utilities for the construction of fermionic interaction operators (represented by the FermionOperator class) and ON vectors (represented by the FermionState class). In the next subsection FermionSpace is discussed, while FermionOperator and FermionState are covered in later sections.

The FermionSpace Class

This class provides a range of tools related to the fermionic Fock space. These tools include the definition of the space itself, and operations on this space related to the generation of fermionic operators and for specification of occupations. The FermionSpace class represents a fermionic Hilbert space, maintaining consistent indexing of spin-orbitals. It incorporates the fermionic anti-commutation algebra in the logic of the class, allowing for the generation of various fermionic quantum operators acting on a given Hilbert space.

Various fermionic quantum operators can be constructed from the methods in this class. These include the number operators, spin operators, one- and two-body operators in various forms, and the anti-Hermitian operators needed for unitary coupled cluster ansatz generation. For example, the UCC terms for minimal basis H₂ (with 4 spin-orbitals) are defined as follows:

from inquanto.spaces import FermionSpace
from inquanto.states import FermionState

# create a FermionSpace object for, say, H2
space = FermionSpace(4)
state = FermionState([1, 1, 0, 0])

# with the state and space objects we can construct operators using native functions, for example:
singles = space.construct_single_ucc_operators(state)
doubles = space.construct_double_ucc_operators(state)
uccsd_excitations = singles + doubles

print("number of uccsd excitations for minimal basis h2:", len(uccsd_excitations))

number of uccsd excitations for minimal basis h2: 3

In this case, the occupations specified by state are needed to construct the UCC operators, though this is not always required depending on the operators being generated. The operators are returned as FermionOperator objects (FermionOperator and FermionState are discussed in more detail later). Several types of operators can be generated by the various methods of FermionSpace, including various forms of number operators, excitation operators and symmetry operators. An exhaustive list of these operators can be found in the API reference.

Once the space and occupations are defined, information relevant to the ON vector can be displayed using the print_state() method:

from inquanto.spaces import FermionSpace
from inquanto.states import FermionState

# create a FermionSpace object for, say, H2
space = FermionSpace(4)
state = FermionState([1, 1, 0, 0])

print("fermionic fock state:")
space.print_state(state)

fermionic fock state:
 0 0a         :  1    
 1 0b         :  1    
 2 1a         :  0    
 3 1b         :  0    

The left-most column refers to spin-orbital indexes, the next column shows the spatial orbital indexes with the alpha (a) and beta (b) spin labels, while the right-most column shows the occupation of each spin-orbital.

Typically, operators generated using the FermionSpace methods will preserve particle number and spin conservation symmetries. In addition to these, FermionSpace can also optionally handle point group symmetry information. By passing in point group symmetries when instantiating the FermionSpace class object, we can generate only those excitations that are non-zero by symmetry. Below, the single excitations of H₂ are symmetry forbidden, leaving only one double.

from inquanto.spaces import FermionSpace
from inquanto.states import FermionState
from inquanto.symmetry import PointGroup

# create a FermionSpace object for, say, H2
space = FermionSpace(4)
state = FermionState([1, 1, 0, 0])

# we can also pass symmetry information to remove redundant excitations
space = FermionSpace(
    4, point_group=PointGroup("D2h"), orb_irreps=["Ag", "Ag", "B1u", "B1u"]
)
singles = space.construct_single_ucc_operators(state)
doubles = space.construct_double_ucc_operators(state)
uccsd_excitations = singles + doubles

print("number of allowed uccsd excitations for minimal basis h2:", len(uccsd_excitations))

number of allowed uccsd excitations for minimal basis h2: 1

This potentially reduces the size of the problem; for example, the number of parameters in variational algorithms like VQE. When a FermionSpace object is generated from a driver, it will typically include point group information if the driver utilized this information in performing precursor classical computation. The FermionSpace classes can be used to generate explicit operator representations of the symmetries of the system with the symmetry_operators_z2() method - usage of this is detailed in the symmetry section.

Fermion Operators & States

In second quantization, operators acting on fermionic states are represented by linear combinations and tensor products of creation and annihilation operators. These obey the anticommutation relations

(71)\[ \begin{align}\begin{aligned}\{\hat{f}^\dagger_P, \hat{f}^\dagger_{P'}\} &= 0\\\{\hat{f}_P, \hat{f}_{P'}\} &= 0\\\{\hat{f}^\dagger_P, \hat{f}_{P'}\} &= \delta_{P, P'}\end{aligned}\end{align} \]

Where \(\{\hat{A}, \hat{B}\} = \hat{A}\hat{B} + \hat{B}\hat{A}\). Mathematically, these operators couple ON vectors belonging to different subspaces \(\mathcal{F}(Q, N)\). Thus, a given arrangement of occupations can be used to build an ON vector by applying a product of creation operators to the vacuum state, such that creation operators are only applied to spin-orbitals or modes that should be unoccupied

(72)\[|\boldsymbol{\eta}\rangle = \prod_{P=0}^{Q-1} (\hat{f}^\dagger_P)^{\eta_P} |0_{0}, ... , 0_{P}, ... , 0_{Q-1}\rangle\]

In InQuanto, vectors in \(\mathcal{F}(Q, N)\) space are represented by FermionState objects. Individual basis states without coefficients are represented by FermionStateString objects. A FermionStateString functions as a dictionary mapping spin-orbital indices represented as integers, to occupation numbers (also represented as integers, of values 0 or 1). For example, for an ON vector of four spin-orbitals, with spin-orbitals indexed 0 and 1 occupied, and spin-orbitals indexed 3 and 4 unoccupied:

from inquanto.states import FermionStateString
on_vector = FermionStateString({0:1, 1:1, 3:0, 4:0})
print(on_vector)

{0: 1, 1: 1, 3: 0, 4: 0}

A simple list of occupation numbers may also be provided to the constructor, as a convenient alternative. In this case, a default register of fermionic modes will be assumed, indexed from \([0,N)\), where \(N\) is the length of the provided list (i.e. the number of spin-orbitals)

from inquanto.states import FermionStateString
on_vector = FermionStateString([1,1,0,0])
print(on_vector)

{0: 1, 1: 1, 2: 0, 3: 0}

Note

The use of a dictionary with integer keys – as opposed to a simple list or tuple – allows for the representation of spin-orbitals with discontinuous indices. This may be useful when representing states which are elements of different Hilbert spaces. It also allows for a uniform interface with classes which represent states wherein the modes are represented by more complex objects, as in the case of QubitState below.

The FermionState then represents the state as a dictionary, with the FermionStateString objects (occupation configurations) as keys, and numerical quantities (configuration coefficients) as values. Typically, this will be constructed with a dict giving such a mapping:

from inquanto.states import FermionState, FermionStateString

on_vector = FermionStateString({0:1, 1:1, 2:0, 3:0})
state = FermionState({on_vector: 1.0})
print(state)

(1.0, {0: 1, 1: 1, 2: 0, 3: 0})

If a simple FermionState representing a single basis state with unit coefficient is required, a convenient alternative for construction is also provided. A list of integer occupation values may be directly passed to the constructor, following the logic of the FermionStateString above:

state = FermionState([1, 1, 0, 0])
print(state)

(1.0, {0: 1, 1: 1, 2: 0, 3: 0})

A FermionState object can also be generated as a method of FermionSpace, by passing the occupation list as an argument. This allows for consistent spin-orbital indexing determined by the FermionSpace instance.

from inquanto.spaces import FermionSpace
from inquanto.states import FermionState

space = FermionSpace(4)
state = space.generate_occupation_state_from_list([1, 1, 0, 0])

FermionState objects may be operated on with standard linear algebraic operations to return new FermionState objects, and may be iterated over.

new_state = 2.0 * state + FermionState([0,0,1,1])
for x in new_state.split():
    print(x)

(2.0, {0: 1, 1: 1, 2: 0, 3: 0})
(1.0, {0: 0, 1: 0, 2: 1, 3: 1})

In addition to representing fermionic states, InQuanto has capabilities for representing operators acting within fermionic Hilbert spaces. Analogously to FermionStateString, the FermionOperatorString represents a single unweighted string of fermionic creation and annihilation operators. A FermionOperatorString may be created with a tuple of pairs of integers, in which the first integer represents the spin-orbital or mode and the second integer is 1 (0) for creation (annihilation) operators. For example, consider an excitation operator in which an electron is annihilated from spin-orbital 2 and created on spin-orbital 3 ( \(\hat{f}^\dagger_3\hat{f}_2\) ):

from inquanto.operators import FermionOperator, FermionOperatorString

f_op_str = FermionOperatorString(((3, 1), (2, 0)))
print(f_op_str)

F3^ F2 

For convenience, it is also possible to create such an object from string input:

f_op_str = FermionOperatorString.from_string("F3^ F2")
print(f_op_str)

F3^ F2 

Similarly to FermionState, fermionic operators are then represented by FermionOperator. This stores the coefficients of each sequence of creation/annihilation operators as a dictionary, in which the keys are FermionOperatorString objects and the values are their coefficients. It may be created through providing the constructor with such a dictionary. Again, a convenient alternative is provided for the construction of an operator comprised of a single string:

f_op_string = FermionOperatorString(((3, 1), (2, 0)))
f_op1 = FermionOperator({f_op_string: 1.0})
f_op2 = FermionOperator(f_op_string, 1.0)
print(f_op1)
print(f_op2)

(1.0, F3^ F2 )
(1.0, F3^ F2 )

Additional construction methods are also available in the reference documentation for FermionOperator. Standard algebraic manipulation is possible for this class:

f_op_string = FermionOperatorString(((3, 1), (2, 0)))
f_op_string_conj = FermionOperatorString(((2, 1), (3, 0)))
f_op = FermionOperator({f_op_string: 1.0})
f_op_conj = FermionOperator({f_op_string_conj: 1.0})
f_op_sum = f_op + f_op_conj
print(f_op_sum)

(1.0, F3^ F2 ), (1.0, F2^ F3 )

Various other methods for manipulating FermionOperator instances are available in InQuanto. The following example ( fermion_operator) summarizes some of the features of InQuanto’s FermionOperator and FermionState objects and their usage.

from inquanto.operators import FermionOperator
from inquanto.operators import FermionOperatorString

# construct a fermion operator
op1 = FermionOperator({FermionOperatorString([(0, 0)]): 1.0})

# now multiply by its adjoint
op = op1 * op1.dagger()

# does it commute with itself? Yes!
print("commutator of op with itself:", op.commutator(op))
print("does op commute with itself?", op.commutes_with(op))

# instantiate from string
op2 = FermionOperator({FermionOperatorString.from_string("F1 F2^"): 3.5})

# is this normal ordered?
print("is op2 normal ordered?", op2.is_normal_ordered())
op2 = op2.normal_ordered()
print("what about now?", op2.is_normal_ordered())
print("is operator number conserving?", op2.is_two_body_number_conserving())

# sum of operators so far:
op3 = op + op2
print("op3 =", op3)
# remove terms with coefficients of absolute value < 3
print("truncated op3 =", op3.truncated(3.0))

# map this fermion operator to a qubit operator
print("JW mapped op3 =", op3.qubit_encode())

# we can apply the operators kets and bras defined in the fermion space to retrieve a new FermionState
from inquanto.states import FermionState

fock_state = FermionState([1, 1, 0, 0])
print("<HF|op3 = ", op3.apply_bra(fock_state))
print("op3|HF> = ", op3.apply_ket(fock_state))

commutator of op with itself: (0.0, ), (0.0, F0^ F0 )
does op commute with itself? True
is op2 normal ordered? False
what about now? True
is operator number conserving? True
op3 = (1.0, F0  F0^), (-3.5, F2^ F1 )
truncated op3 = (-3.5, F2^ F1 )
JW mapped op3 = (0.5, ), (0.5, Z0), (-0.875j, Y1 X2), (-0.875, X1 X2), (-0.875, Y1 Y2), (0.875j, X1 Y2)
<HF|op3 =  (0)
op3|HF> =  (-3.5, {0: 1, 1: 0, 2: 1, 3: 0})

InQuanto also contains a class for dealing with sets of fermionic operators that are not linearly combined – the FermionOperatorList. Each separate FermionOperator in the FermionOperatorList may be associated with an additional scalar. This is particularly useful when describing sequences of exponentiated fermionic operators, for example when considering variational ansatzes.

from sympy import sympify
from inquanto.operators import FermionOperatorList
fol = FermionOperatorList([(sympify("a"),op1),(sympify("b"),op2)])
print(fol)

a         [(1.0, F0 )],
b         [(-3.5, F2^ F1 )]

Note

FermionOperatorList are designed to retain ordering for use in cases where operator ordering matters, such as exponentiation. In contrast, FermionOperator is not naturally ordered.

This class is particularly useful for containing Trotter sequences. The FermionOperator class contains various helper methods for generating such sequences.

op1 = FermionOperator({FermionOperatorString([(0, 0),(1,1)]): 1.0})
op2 = FermionOperator({FermionOperatorString([(2, 0),(3,1)]): 1.0})
op3 = op1 + op2
op_trotterized = op3.trotterize(trotter_number=2)
print(op_trotterized)

0.5       [(1.0, F0  F1^)],
0.5       [(1.0, F2  F3^)],
0.5       [(1.0, F0  F1^)],
0.5       [(1.0, F2  F3^)]

Additional information on the functionality of these classes is provided in the API reference.

Qubit spaces, operators and states

For the analysis of many quantum algorithms, it is useful to work in a representation above the level of quantum circuits decomposed into some set of quantum gate primitives. Similarly to the FermionSpace class, qubit Hilbert spaces are represented in InQuanto using the QubitSpace class.

from inquanto.spaces import QubitSpace
qubit_space = QubitSpace(4)

The QubitSpace object represents the \(2^N\) dimensional Hilbert space of an \(N\)-qubit system. For chemistry purposes, operators and states in a qubit Hilbert space are typically generated from fermionic operators and states, and thus the QubitSpace class does not have as many methods for directly generating operators and states as the fermionic equivalent. However, methods for finding symmetry operators of a given qubit operator are included, as discussed in the symmetry section.

Qubit Operators

Conventionally, operators acting upon states of qubits are typically represented as linear combinations of Pauli strings

(73)\[\hat{O} = \sum_i{h_i} P_i\]

where each Pauli string

(74)\[P_i = \bigotimes^N_{n=0}{p_n}\]

is a tensor product of single qubit Pauli \(p_n \in \{I,X,Y,Z\}\) operators and \(N\) is the number of qubits. Any operator acting on a register of qubits can be decomposed in this form. This representation is particularly useful when considering quantum chemistry problems, where operators acting on states of fermions must be mapped to operators acting on states of qubits prior to simulation. In InQuanto, operators of this form are implemented using the QubitOperator class. In chemistry, instances of this will frequently derive from mapping a FermionOperator, however they may also be constructed manually. Akin to the fermion operators discussed above, these comprise a map between QubitOperatorString Pauli strings and numerical or symbolic coefficients.

Whereas FermionOperatorString represents strings of creation and annihilation operators as a collection of pairs of integers (indexing modes) and \(0\) or \(1\) (annihilation/creation), pytket provides us a more sophisticated approach to identifying qubits. As such, a QubitOperatorString contains a map between pytket qubit objects and Pauli objects, indicating which Pauli operator (of \(\{I,X,Y,Z\}\)) is acting on which qubit. Note that this map is not necessarily exhaustive – while it may contain explicit identity operators, the operator is assumed to also implicitly act with the identity on any qubit not contained within the map.

from inquanto.operators import QubitOperatorString
from pytket import Qubit
from pytket.pauli import Pauli

pauli_string = QubitOperatorString({Qubit(0):Pauli.X, Qubit(1):Pauli.X})
print(pauli_string.map)

{q[0]: <Pauli.X: 1>, q[1]: <Pauli.X: 1>}

Much like FermionOperator, a QubitOperator then stores a linear combination of Pauli strings as a map between these QubitOperatorString objects and numerical or symbolic coefficients.

from inquanto.operators import QubitOperator
pauli_string = QubitOperatorString({Qubit(0):Pauli.X, Qubit(1):Pauli.X})
op = QubitOperator({pauli_string:1.0})
print(op)

(1.0, X0 X1)

As with their fermionic counterparts, a variety of convenient alternative methods to construct the operator are provided, detailed in the API reference for QubitOperator. Basic linear algebra is supported on these objects, along with various useful tools detailed in the API reference for QubitOperator.

op1 = QubitOperator("X0",1.0)
op2 = QubitOperator("Z0",1.0)
op3 = op1 + op2
print('Op1 + Op1:', op3)
print('Op1 conjugated:', op1.dagger())
print('Commutator:', op1.commutator(op2))
print('Sparse matrix form:', op1.to_sparse_matrix(1))
print('Identity:', QubitOperator.identity())

Op1 + Op1: (1.0, X0), (1.0, Z0)
Op1 conjugated: (1.0, X0)
Commutator: (-2j, Y0)
Sparse matrix form: <Compressed Sparse Column sparse matrix of dtype 'complex128'
	with 2 stored elements and shape (2, 2)>
  Coords	Values
  (1, 0)	(1+0j)
  (0, 1)	(1+0j)
Identity: (1.0, )

Note

These examples have used an alternative string construction method for the QubitOperator, which is useful for simple examples such as this.

Several helper methods for determining properties of an operator are also available:

print('Is operator Hermitian?', op3.is_hermitian())
print('Is operator anti-Hermitian?', op3.is_antihermitian())
print('Is operator unitary?', op3.is_unitary())
print('Is operator unit-norm?', op3.is_unit_norm(order=2))
print('Is operator self-inverse?', op3.is_self_inverse())

Is operator Hermitian? True
Is operator anti-Hermitian? False
Is operator unitary? False
Is operator unit-norm? False
Is operator self-inverse? False

As mentioned above, a QubitOperator consists of a mapping between QubitOperatorString objects and their corresponding coefficients in the linear combination representing the operator. The terms may be extracted with the pauli_strings() property, whereas the coefficients may be extracted as a list with the coefficients() property.

print('Operator coefficients:')
print(op3.coefficients)
print('Pauli strings:')
print(op3.pauli_strings)

Operator coefficients:
[1.0, 1.0]
Pauli strings:
[(Xq[0]), (Zq[0])]

Often when considering quantum algorithms, it is useful to describe qubit operators in their symplectic representation. This consists of an \((M\times 2N)\) binary matrix where \(M\) is the number of independent Pauli strings in the operator, and \(N\) is the number of qubits. The leftmost half of this matrix designates whether a Pauli \(X\) is acting on qubit \(n\) in term \(m\), and the rightmost half designates whether a Pauli \(Z\) is acting on the same qubit in the same term. As \(Y = \mathrm{i} XZ\), both leftmost and rightmost entries are \(1\) if a Pauli \(Y\) is present. Note that in this representation, information regarding the coefficients of each term must be stored independently.

print(op1.symplectic_representation())

[[ True False]]

InQuanto additionally provides a class for dealing with sets of Pauli operators that are not linearly combined, akin to the FermionOperatorList – the QubitOperatorList. This provides similar functionality to that of its fermionic counterpart.

from sympy import sympify
from inquanto.operators import QubitOperatorList
qol = QubitOperatorList([(sympify("a"),op1),(sympify("b"),op2)])
print(qol)

a         [(1.0, X0)],
b         [(1.0, Z0)]

Trotterization functionality is also available for the QubitOperator.

op1 = QubitOperator("X0",1.0)
op2 = QubitOperator("Z0",1.0)
op3 = op1 + op2
op_trotterized = op3.trotterize(trotter_number=2)
print(op_trotterized)

0.5       [(1.0, X0)],
0.5       [(1.0, Z0)],
0.5       [(1.0, X0)],
0.5       [(1.0, Z0)]

The API reference details these methods, and provides a full breakdown of the other functionality available for the QubitOperatorList classes.

Qubit States & Expectation Values

A register of \(N\) qubits corresponds to a \(\mathcal{C}^{2^{N}}\) Hilbert space. As such, it can be represented with a \(2^N\) dimensional vector of complex numbers. Clearly, generating such a vector is not scalable on a classical computer (otherwise we wouldn’t need quantum computers), but this approach can be practical for small \(N\). For some tasks (for instance, if we are simulating the action of a known Clifford operator), we can guarantee that a given state will have polynomial support. In this case, we can efficiently store the state in a sparse state vector. InQuanto provides an alternative to an explicit sparse state vector representation of states in the form of the QubitState class, instances of which consist of linear combinations of QubitStateString objects.

from inquanto.states import QubitState
qubit_state = QubitState([1,1,0,0],1.)
print(qubit_state)

(1.0, {q[0]: 1, q[1]: 1, q[2]: 0, q[3]: 0})

These can be converted to state vector representations:

state_vector = qubit_state.to_ndarray()
print(state_vector)

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

Most importantly for the purposes of analyzing quantum algorithms, the QubitState representation allows for performing linear algebra with other QubitState and QubitOperator objects.

overlap = qubit_state.vdot(qubit_state)
expectation_value = op1.state_expectation(qubit_state)
print('Overlap:', overlap)
print('Expectation value:', expectation_value)

Overlap: 1.0
Expectation value: 0

QubitState objects are implicitly sparse, in contrast to full dense statevector representations. Performing calculations in this way allows for the analysis of qubit states and operators without the need to either generate full circuit representations, or the expensive generation of full \(2^{2N}\) matrix representations of operators.

Fermion-to-Qubit Mapping

In quantum chemistry, we are most often concerned with the properties of electrons. As electrons are fermions, second-quantized fermionic creation and annihilation operators obey the fermionic anticommutation relations:

(75)\[\{\hat{a}^\dagger_i , \hat{a}^\dagger_j\} = 0, \{\hat{a}_i , \hat{a}_j\} = 0, \{\hat{a}^\dagger_i , \hat{a}_j\} = \delta_{i,j}\]

This implicitly restricts the occupation number of a given fermionic mode to \(\{0,1\}\). Conversely, qubits can be considered to be paraparticles. Like fermionic modes, each qubit is a two-level system. However, the above fermionic anticommutation relations are not obeyed:

(76)\[\begin{split}\mathinner{|{1}\rangle} \mathinner{\langle{0}|}_i = X_i - \mathrm{i}Y_i, \\ \mathinner{|{0}\rangle} \mathinner{\langle{1}|}_i = X_i + \mathrm{i}Y_i \\ \{X_i - \mathrm{i}Y_i , X_j + \mathrm{i}Y_j\} \neq \delta_{i,j}\end{split}\]

As such, in order to use a system of qubits to simulate a system of fermions, the fermionic anticommutation relations must be encoded. An encoding scheme consists of a linear map between states and operators in the fermionic Hilbert space and the states and operators in the qubit Hilbert space. For quantum chemistry purposes, mapping both states and operators is important - for example, in a canonical VQE calculation, both the Hamiltonian operator and the reference state must be mapped to the qubit space.

In InQuanto, fermion-to-qubit mappings are stored in the inquanto.mappings module. Four mappings are included in the current version:

• Jordan-Wigner transformation

• Bravyi-Kitaev mapping

• Parity mapping

• "Paraparticular" mapping – this does not encode fermionic statistics

In order to map an InQuanto FermionOperator to a QubitOperator, the operator_map() method can be used:

from inquanto.operators import FermionOperator,FermionOperatorString
from inquanto.mappings import QubitMappingJordanWigner
fermion_operator = FermionOperator(FermionOperatorString([(1, 0)]), 0.5)
qubit_operator = QubitMappingJordanWigner.operator_map(fermion_operator)
print(qubit_operator)

(0.25, Z0 X1), (0.25j, Z0 Y1)

and similarly, to map FermionState objects to QubitState objects, the state_map() method is used:

from inquanto.states import FermionState
fermion_state = FermionState([0,0,1,1])
qubit_state = QubitMappingJordanWigner.state_map(fermion_state)

Note that both of these methods include an optional argument specifying the register of tket Qubit objects that comprise the target qubit space. Without this argument provided (as above), the mapping will if possible infer that a minimally sized register is to be used, indexed incrementally from 0. For some mappings (for example, the Bravyi-Kitaev mapping), this information cannot be inferred. If these mappings are used, it is necessary to specify the qubit register.

from pytket import Qubit
from inquanto.mappings import QubitMappingBravyiKitaev
qubit_operator_bk = QubitMappingBravyiKitaev.operator_map(fermion_operator,[Qubit(i) for i in range(8)])
print(qubit_operator_bk)

(0.25, Z0 X1 X3 X7), (0.25j, Y1 X3 X7)

Finally, when performing state vector simulations one may wish to represent fermionic and qubit states in the form of full dense or sparse vectors of complex numbers. The inquanto.mappings classes support this, and will return in the same representation as the input.

from numpy import array
from scipy.sparse import csc_matrix
dense_fermionic_state = array([[1],[0],[0],[0]])
sparse_fermionic_state = csc_matrix(dense_fermionic_state)
dense_qubit_state = QubitMappingJordanWigner.state_map(dense_fermionic_state,qubits=[Qubit(i) for i in range(4)])
sparse_qubit_state = QubitMappingJordanWigner.state_map(sparse_fermionic_state,qubits=[Qubit(i) for i in range(4)])
print(dense_qubit_state)
print(sparse_qubit_state)

[[1]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]]
<Compressed Sparse Column sparse matrix of dtype 'int64'
	with 1 stored elements and shape (16, 1)>
  Coords	Values
  (0, 0)	1

Custom mappings

For advanced usage, it is possible to define custom mapping schemes. The QubitMapping class is a superclass containing the logic of the operator_map() and state_map() methods. Custom mappings can be designed by subclassing this and implementing the update_set() , parity_set() , rho_set() and state_map_matrix() methods. The first three of these methods return sets of qubits according to the formalism of Seeley, Richard and Love [38]. The state_map_matrix() method returns a matrix which transforms a binary column vector representation of the binary index corresponding to a given basis state in the fermionic space, to a similar vector in the qubit space. In the flip/update/parity/rho set formalism, fermionic creation and annihilation operators are mapped according to the following relation:

(77)\[ \begin{align}\begin{aligned}&a^\dagger_i \rightarrow \frac{1}{2} (X_{U(i)} X_i Z_{P(i)} - \mathrm{i}X_{U(i)} Y_i Z_{\rho(i)})\\&a_i \rightarrow \frac{1}{2} (X_{U(i)} X_i Z_{P(i)} + \mathrm{i}X_{U(i)} Y_i Z_{\rho(i)})\end{aligned}\end{align} \]

where \(U(i)\) is the update set, \(P(i)\) is the parity set, \(\rho(i)\) is the rho set for orbital \(i\). As an example, for the Jordan-Wigner mapping:

(78)\[ \begin{align}\begin{aligned}&a^\dagger_i \rightarrow \frac{1}{2} \left(X_i - \mathrm{i} Y_i \right) \otimes Z^{\otimes i}\\&a_i \rightarrow \frac{1}{2} \left(X_i + \mathrm{i} Y_i \right) \otimes Z^{\otimes i}\end{aligned}\end{align} \]

Comparing this against (77) we see the update and flip set are always empty, whereas the parity set consists of all qubits with index less than \(i\). state_map_matrix() for the Jordan-Wigner transformation will always return the identity matrix, since state \(\ket{n}\) in the fermionic space is always mapped to state \(\ket{n}\) in the qubit space.

Integral Operators

We have described above how InQuanto stores and handles the fermionic Hamiltonian via the FermionOperator class. The FermionOperator class stores those operators and numeric integral values as items of a dictionary. Such storage is convenient for manipulating and accessing individual terms of the Hamiltonian, but is not tailored well to linear algebra operations, such as integral transformation.

The InQuanto operators module provides integral operator classes for storing and manipulating chemistry integrals (\(h_{ij}\) and \(h_{ijkl}\), typically molecular orbital integrals) as algebraic objects. The standard integral operators are the ChemistryRestrictedIntegralOperator and ChemistryUnrestrictedIntegralOperator classes, for spin-restricted and spin-unrestricted formalisms respectively. In the restricted case, one-body integrals are stored in a 2-dimensional numpy array with shape (n, n), and two-body integrals in a 4-dimensional array with shape (n, n, n, n), where n is the number of spatial orbitals. In the unrestricted case, the \(\uparrow\) and \(\downarrow\) spin channels have independent integrals. Hence, two 2-dimensional arrays store the one-body integrals, one for each spin channel, and four 4-dimensional arrays store all spin-configurations of the two-body integrals (\(\upuparrows\upuparrows, \downdownarrows\downdownarrows, \upuparrows\downdownarrows\) and \(\downdownarrows\upuparrows\)).

The inquanto-pyscf extension is the primary tool for generating integral operators for a chemical system, though integral operators may also be instantiated directly with numpy arrays (see the ChemistryRestrictedIntegralOperator and ChemistryUnrestrictedIntegralOperator constructors). Below, we use the express module to load in a pre-computed integral operator for LiH:

import pandas as pd
from inquanto.express import load_h5

pd.options.display.max_rows = 15
lih_sto3g = load_h5('lih_sto3g.h5', as_tuple = True)
integral_operator = lih_sto3g.hamiltonian_operator
print(integral_operator.df())

      Coefficients                Terms
0         1.050650                     
1        -4.746695              F0^ F0 
2         0.109103              F0^ F2 
3         0.168024              F0^ F4 
4        -0.026783             F0^ F10 
...            ...                  ...
1496      0.009754    F11^ F10^ F8  F9 
1497      0.009754    F11^ F10^ F8  F9 
1498      0.228088  F11^ F10^ F10  F11 
1499      0.228088  F11^ F10^ F10  F11 
1500     -0.935480            F11^ F11 

[1501 rows x 2 columns]

where the df() method produces a pandas dataframe of all integrals and their corresponding fermion operator terms. Integral operators may be converted directly into a qubit Hamiltonian with the qubit_encode() method, which uses Jordan-Wigner mapping by default:

qubit_hamiltonian = integral_operator.qubit_encode()
print("LiH STO-3G JW qubit hamiltonian:\n", qubit_hamiltonian.df())

LiH STO-3G JW qubit hamiltonian:
      Coefficient                                   Term  \
0      -4.107196                                          
1      -0.396888                                    Z11   
2      -0.396888                                    Z10   
3       0.114044                                Z10 Z11   
4      -0.228984                                     Z9   
..           ...                                    ...   
626    -0.011161  Z0 X1 Z2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 X11   
627     0.029058                            Z0 Y1 Z2 Y3   
628     0.034437                      Z0 Y1 Z2 Z3 Z4 Y5   
629    -0.011161  Z0 Y1 Z2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Y11   
630     0.414556                                  Z0 Z1   

            Coefficient Type  
0    <class 'numpy.float64'>  
1    <class 'numpy.float64'>  
2    <class 'numpy.float64'>  
3    <class 'numpy.float64'>  
4    <class 'numpy.float64'>  
..                       ...  
626  <class 'numpy.float64'>  
627  <class 'numpy.float64'>  
628  <class 'numpy.float64'>  
629  <class 'numpy.float64'>  
630  <class 'numpy.float64'>  

[631 rows x 3 columns]

Unitary transformations may be applied to the integrals with the rotate() method, which takes a unitary matrix and transforms in-place:

from scipy.stats import ortho_group
u = ortho_group.rvs(lih_sto3g.n_orbital)  # Random, real unitary matrix
integral_operator.rotate(u)
print(integral_operator.df())

      Coefficients                Terms
0         1.050650                     
1        -1.151561              F0^ F0 
2         0.189985              F0^ F2 
3        -0.051761              F0^ F4 
4         0.032263              F0^ F6 
...            ...                  ...
4460      0.014764   F11^ F10^ F9  F10 
4461      0.014764   F11^ F10^ F9  F10 
4462      0.268396  F11^ F10^ F10  F11 
4463      0.268396  F11^ F10^ F10  F11 
4464     -1.324213            F11^ F11 

[4465 rows x 2 columns]

Integral operators support other utility methods including approx_equal(), for comparing Hamiltonians, items(), which iterates over terms, and to_FermionOperator(), for converting between operator and integral focused objects.

Alongside integral operators, the operators module also provides a set of classes for managing reduced density matrices (RDMs). These may be used in combination with integral operators to compute useful properties of the Hamiltonian. For example, with the UnrestrictedOneBodyRDM we may compute the total mean-field energy and effective potential matrices:

from inquanto.operators import UnrestrictedOneBodyRDM
import numpy as np

integral_operator = load_h5('h3_sto3g_m2_u.h5', as_tuple = True).hamiltonian_operator
rdm1 = UnrestrictedOneBodyRDM(rdm1_aa=np.diag([1, 1, 0]), rdm1_bb=np.diag([1, 0, 0]))

print("Mean-field energy:\n", integral_operator.energy(rdm1))
print("\nEffective potential a:\n", integral_operator.effective_potential(rdm1)[0])
print("\nEffective potential b:\n", integral_operator.effective_potential(rdm1)[1])

Mean-field energy:
 -1.5140974066187696

Effective potential a:
 [[ 0.996 -0.     0.226]
 [-0.     0.924  0.   ]
 [ 0.226  0.     1.632]]

Effective potential b:
 [[ 1.163 -0.     0.069]
 [-0.     1.541 -0.   ]
 [ 0.069  0.     1.712]]

Above, RDM is initialized with numpy arrays in the same basis as the integral operator. The energy() method may also be provided with a 2-RDM (see RestrictedTwoBodyRDM and UnrestrictedTwoBodyRDM) to calculate the total, non-mean-field, energy.

In addition to the basic integral operator classes, InQuanto also includes compact integral operators: ChemistryRestrictedIntegralOperatorCompact and ChemistryUnrestrictedIntegralOperatorCompact, which exploit symmetries in the two-body integrals to reduce classical memory requirements. Compact integral operator classes support the same operations as discussed above, and are most naturally instantiated using the inquanto-pyscf extension.

Orbital Transformation and Optimization

Slater determinants, represented by vectors in Fock space, are functions of molecular orbitals. Typically the molecular orbitals are found by solving some other problem, such as the Hartree-Fock equations, but in some cases orbitals are chosen to be some other set of functions. For instance, they could be localized or optimized in some other way. InQuanto contains tools to help with procedures of this type, and to facilitate the transfer of these methods to quantum algorithms.

The OrbitalTransformer class contains several methods for common practices in molecular orbital manipulation. For instance, to Gram-Schmidt orthogonalize a set of molecular orbitals in an orthogonal atomic orbital basis,

from inquanto.operators import OrbitalTransformer
import numpy
orbitals = numpy.array([[1, 2], [3, 4]])
ot = OrbitalTransformer()
ot.gram_schmidt(v=orbitals, overlap=None)

array([[-0.316, -0.949],
       [-0.949,  0.316]])

or, equivalently, in a non-orthogonal atomic orbital basis,

s = numpy.array([[1.0, 0.66314574], [0.66314574, 1.0]])
ot.gram_schmidt(v=orbitals, overlap=s)

array([[-0.267, -1.309],
       [-0.802,  1.068]])

One can achieve the same goal of orthonormalization using the orthonormalize() method, which finds the closest orthonormal set by the symmetric transformation.

ot.orthonormalize(v=orbitals, overlap=s)

array([[-0.887,  0.999],
       [ 1.336,  0.001]])

The OrbitalTransformer object also defines a method for computing the unitary which relates two sets of molecular orbitals. For example, to find the unitary relating the MO coefficients in the matrix \(X\) with those in matrix \(C\),

X = numpy.array([[1, 2], [3, 4]])
C = numpy.array([[2, 1], [4, 3]])
ot = OrbitalTransformer()
my_unitary = ot.compute_unitary(v_init=X, v_final=C)
print(my_unitary)

[[0. 1.]
 [1. 0.]]

Similarly, to transform orbitals X to orbitals C, if the unitary is known,

new_C = ot.transform(v=X, tu=my_unitary)
print(new_C)

[[2. 1.]
 [4. 3.]]

The majority of the time, the unitary is not known, and is the result of some optimization process. The OrbitalOptimizer class in InQuanto is constructed with a black-box function which finds the rotational unitary, given some variational criteria. For instance, we can perform a localization which depends on some function localize:

from inquanto.operators import OrbitalOptimizer
oo = OrbitalOptimizer(
     v_init=initial_orbitals,
     occ=[2, 2, 0, 0],
     split_rotation=False,
     functional=localize,
     minimizer=MinimizerScipy(),
     reduce_free_parameters=True
)
final_orbitals, minimising_unitary, final_value = oo.optimize()

The OrbitalOptimizer object will try to retain orbital symmetries if point_group and orb_irreps are both passed into the constructor. The function passed to the functional argument must be a function of a 2D array.

Note

The OrbitalOptimizer is also a callable object, but the callable execution of the optimization returns only the optimized orbitals. This is for compatibility with the transf functionality in some extensions.

After the optimization, a report can be generated with the generate_report() method.

Double Factorization

An important restriction in near-term chemistry applications of quantum computing is the size of the Hamiltonian operator, particularly the two-body interaction term:

(79)\[\hat H_2 = \frac{1}{2}\sum_{pqrs} h_{pqrs} a_p^\dagger a_r^\dagger a_s a_q,\]

which has \(\mathcal{O}(N^4)\) terms where \(N\) is the number of orbitals. Double factorization is a two-step, tensor decomposition of the two-body integrals \(h_{pqrs}\), which provides a systematic approach for truncating terms in two-body Hamiltonians. Below, we provide an introduction to the essential equations, and an example of using this decomposition strategy with InQuanto.

The first decomposition is given by:

(80)\[h_{pqrs} = \sum_t^{N_\gamma} V_{pq}^t \gamma^t V_{rs}^t,\]

where \(N_\gamma = N^2\) for an exact decomposition. And the second:

(81)\[V_{pq}^t = \sum_u^{N_\lambda^t} U_{pu}^t \lambda_u^t U_{qu}^t.\]

for each \(t\), where \(N_\lambda^t = N\) for an exact decomposition.

In InQuanto, the first factorization can be done either by eigenvalue decomposition or a pivoted, incomplete Cholesky decomposition [39, 40]. The second factorization step is always an eigenvalue decomposition, since the fully factorized operator must be expressed in terms of unitary matrices to enable circuit-level rotations as shown below.

Each of the decompositions may be truncated to better control the size and scaling of number of terms in the Hamiltonian. The magnitudes of the factors \(\gamma^t\) and \(\lambda_u^t\) determine the importance of their corresponding terms in the decompositions and so terms can be systematically discarded based on their values.

The sums are truncated by discarding terms in ascending order of the factor magnitudes starting from the smallest such that the sum of the discarded factors does not exceed some threshold [41]:

(82)\[\sum_{t=N_\gamma +1}^{N^2} \left|\gamma^t\right| < \varepsilon_1,\,\,\,\,\sum_{u=N_\lambda^t +1}^{N} \left|\lambda^t_u\right| < \varepsilon_2.\]

In practice, the first decomposition can be truncated to \(N_\gamma\sim \mathcal{O}(N)\), and the second decomposition to \(N_\lambda^t < N\), helping to control the scaling of the number of terms in the two-body Hamiltonian.

Inserting the double factorized integrals into the two-body Hamiltonian (79), after some work we arrive at the diagonalized expression [41]:

(83)\[\hat{H}_2 = \hat{S} + \frac{1}{2}\sum_t^{N_\gamma} \gamma^t \hat{R}(\textbf{U}_t) \left( \sum_u^{N_\lambda^t} \lambda_u^t a_u^\dagger a_u \right)^2 \hat{R}(\textbf{U}_t)^\dagger.\]

Here, \(\hat{S}\) is a one-body offset term which comes from rearranging the fermion operators:

(84)\[\hat{S} = -\frac{1}{2}\sum_{pq} \left( \sum_r h_{prrq}\right) a_p^\dagger a_q,\]

and \(\hat{R}(\textbf{U}_t)\) are Fock-space basis rotation operators, corresponding to the single particle basis rotation matrices \(\textbf{U}_t\) (the eigenvector matrices of the second decomposition (81)). These operators are given by the Thouless theorem [42]:

(85)\[\hat{R}(\textbf{A}) = \exp\left[ \sum_{ij} \left[ \ln \textbf{A}\right]_{ij} a_i^\dagger a_j \right],\]

which are implemented at the circuit level with a Givens QR decomposition [43]. See below for an example, and here for more details.

We may similarly diagonalize the one-body integrals and \(\hat{S}\) operator:

(86)\[\hat{H}_1' = \hat{H}_1 + \hat{S} = \sum_{pq} h_{pq}' a_p^\dagger a_q = \hat{R}(\textbf{W})\left( \sum_r \omega_r a_r^\dagger a_r \right)\hat{R}(\textbf{W})^\dagger\]

where all one-body-like terms have been consolidated into an effective one-body Hamiltonian \(H_1'\). The full Hamiltonian can then be written in this diagonal form as:

(87)\[\begin{split}\hat{H} =\, &\hat{H}_0 + \hat{R}(\textbf{W})\left( \sum_r \omega_r a_r^\dagger a_r \right)\hat{R}(\textbf{W})^\dagger \\ &+ \frac{1}{2}\sum_t^{N_\gamma} \gamma^t \hat{R}(\textbf{U}_t) \left( \sum_u^{N_\lambda^t} \lambda_u^t a_u^\dagger a_u \right)^2 \hat{R}(\textbf{U}_t)^\dagger.\end{split}\]

In InQuanto, a ChemistryRestrictedIntegralOperator object may be transformed into this representation, a DoubleFactorizedHamiltonian object, using the double_factorize() method. We demonstrate this below for H₂O in the STO-3G basis:

import numpy as np
from inquanto.express import load_h5
from inquanto.operators import ChemistryRestrictedIntegralOperator

ham = load_h5("h2o_sto3g.h5", as_tuple=True).hamiltonian_operator
df_ham = ham.double_factorize(
    tol1=1e-3,
)

gammas = df_ham.two_body.outer_params
lambdas = df_ham.two_body.inner_params
N_gamma = len(gammas)
N_lambda = [len(l) for l in lambdas]

df_ham.n_orb, N_gamma, N_lambda

(7, 23, [7, 7, 4, 7, 2, 7, 4, 7, 2, 2, 7, 4, 4, 7, 2, 7, 4, 7, 2, 7, 2, 4, 7])

where we show the extent of the truncation (recall equation (82)). Default behavior is to consolidate and diagonalize all one-body terms, so df_ham represents a Hamiltonian of the form (87).

Note

double_factorize() performs the decomposition on the molecular orbital (MO) integrals. In classical literature, similar methods have been used to reduce the memory storage requirements of two-body atomic orbital (AO) integrals [44]. In InQuanto, the purpose of this approach is to truncate the Hamiltonian for quantum simulation, so it should not be expected to reduce classical memory requirements.

The items() method returns an iterator which generates the FermionOperator terms (the number operator sums in brackets of (87)) alongside their corresponding rotation matrices (\(\textbf{W}\) and \(\textbf{U}_t\)). For example, the one-body and the first two-body terms are:

terms = list(df_ham.items())
print("One body operator: ", terms[1][0])
print("One body rotation: ", terms[1][1])
print("\nFirst two body operator: ", terms[2][0])
print("First two body rotation: ", terms[2][1])

One body operator:  (-35.15627858132068, F0^ F0 ), (-35.15627858132068, F1^ F1 ), (-8.99140515493988, F2^ F2 ), (-8.99140515493988, F3^ F3 ), (-8.30673410272335, F4^ F4 ), (-8.30673410272335, F5^ F5 ), (-8.057439359193436, F6^ F6 ), (-8.057439359193436, F7^ F7 ), (-8.028943982634141, F8^ F8 ), (-8.028943982634141, F9^ F9 ), (-4.621602713403972, F10^ F10 ), (-4.621602713403972, F11^ F11 ), (-4.593693566394927, F12^ F12 ), (-4.593693566394927, F13^ F13 )
One body rotation:  [[-0.999 -0.031 -0.    -0.011  0.     0.     0.006]
 [ 0.028 -0.897 -0.    -0.22   0.    -0.    -0.382]
 [-0.     0.    -0.775  0.     0.     0.632 -0.   ]
 [ 0.011  0.046 -0.    -0.908  0.     0.     0.416]
 [ 0.    -0.    -0.    -0.    -1.    -0.     0.   ]
 [ 0.014 -0.438 -0.     0.355 -0.     0.     0.825]
 [-0.     0.    -0.632  0.     0.    -0.775  0.   ]]

First two body operator:  (2.046238405452405, F0^ F0  F0^ F0 ), (2.046238405452405, F0^ F0  F1^ F1 ), (0.8165371440369917, F0^ F0  F2^ F2 ), (0.8165371440369917, F0^ F0  F3^ F3 ), (0.8055041237052253, F0^ F0  F4^ F4 ), (0.8055041237052253, F0^ F0  F5^ F5 ), (0.767389123035518, F0^ F0  F6^ F6 ), (0.767389123035518, F0^ F0  F7^ F7 ), (0.7479034551648777, F0^ F0  F8^ F8 ), (0.7479034551648777, F0^ F0  F9^ F9 ), (0.41264550253745813, F0^ F0  F10^ F10 ), (0.41264550253745813, F0^ F0  F11^ F11 ), (0.4035573895111223, F0^ F0  F12^ F12 ), (0.4035573895111223, F0^ F0  F13^ F13 ), (2.046238405452405, F1^ F1  F0^ F0 ), (2.046238405452405, F1^ F1  F1^ F1 ), (0.8165371440369917, F1^ F1  F2^ F2 ), (0.8165371440369917, F1^ F1  F3^ F3 ), (0.8055041237052253, F1^ F1  F4^ F4 ), (0.8055041237052253, F1^ F1  F5^ F5 ), (0.767389123035518, F1^ F1  F6^ F6 ), (0.767389123035518, F1^ F1  F7^ F7 ), (0.7479034551648777, F1^ F1  F8^ F8 ), (0.7479034551648777, F1^ F1  F9^ F9 ), (0.41264550253745813, F1^ F1  F10^ F10 ), (0.41264550253745813, F1^ F1  F11^ F11 ), (0.4035573895111223, F1^ F1  F12^ F12 ), (0.4035573895111223, F1^ F1  F13^ F13 ), (0.8165371440369917, F2^ F2  F0^ F0 ), (0.8165371440369917, F2^ F2  F1^ F1 ), (0.32583344434133926, F2^ F2  F2^ F2 ), (0.32583344434133926, F2^ F2  F3^ F3 ), (0.32143079463649665, F2^ F2  F4^ F4 ), (0.32143079463649665, F2^ F2  F5^ F5 ), (0.30622126982800785, F2^ F2  F6^ F6 ), (0.30622126982800785, F2^ F2  F7^ F7 ), (0.2984456501590826, F2^ F2  F8^ F8 ), (0.2984456501590826, F2^ F2  F9^ F9 ), (0.16466330572421778, F2^ F2  F10^ F10 ), (0.16466330572421778, F2^ F2  F11^ F11 ), (0.16103675769568102, F2^ F2  F12^ F12 ), (0.16103675769568102, F2^ F2  F13^ F13 ), (0.8165371440369917, F3^ F3  F0^ F0 ), (0.8165371440369917, F3^ F3  F1^ F1 ), (0.32583344434133926, F3^ F3  F2^ F2 ), (0.32583344434133926, F3^ F3  F3^ F3 ), (0.32143079463649665, F3^ F3  F4^ F4 ), (0.32143079463649665, F3^ F3  F5^ F5 ), (0.30622126982800785, F3^ F3  F6^ F6 ), (0.30622126982800785, F3^ F3  F7^ F7 ), (0.2984456501590826, F3^ F3  F8^ F8 ), (0.2984456501590826, F3^ F3  F9^ F9 ), (0.16466330572421778, F3^ F3  F10^ F10 ), (0.16466330572421778, F3^ F3  F11^ F11 ), (0.16103675769568102, F3^ F3  F12^ F12 ), (0.16103675769568102, F3^ F3  F13^ F13 ), (0.8055041237052253, F4^ F4  F0^ F0 ), (0.8055041237052253, F4^ F4  F1^ F1 ), (0.32143079463649665, F4^ F4  F2^ F2 ), (0.32143079463649665, F4^ F4  F3^ F3 ), (0.3170876333750909, F4^ F4  F4^ F4 ), (0.3170876333750909, F4^ F4  F5^ F5 ), (0.30208361911523307, F4^ F4  F6^ F6 ), (0.30208361911523307, F4^ F4  F7^ F7 ), (0.2944130633377986, F4^ F4  F8^ F8 ), (0.2944130633377986, F4^ F4  F9^ F9 ), (0.16243838109804687, F4^ F4  F10^ F10 ), (0.16243838109804687, F4^ F4  F11^ F11 ), (0.15886083485519145, F4^ F4  F12^ F12 ), (0.15886083485519145, F4^ F4  F13^ F13 ), (0.8055041237052253, F5^ F5  F0^ F0 ), (0.8055041237052253, F5^ F5  F1^ F1 ), (0.32143079463649665, F5^ F5  F2^ F2 ), (0.32143079463649665, F5^ F5  F3^ F3 ), (0.3170876333750909, F5^ F5  F4^ F4 ), (0.3170876333750909, F5^ F5  F5^ F5 ), (0.30208361911523307, F5^ F5  F6^ F6 ), (0.30208361911523307, F5^ F5  F7^ F7 ), (0.2944130633377986, F5^ F5  F8^ F8 ), (0.2944130633377986, F5^ F5  F9^ F9 ), (0.16243838109804687, F5^ F5  F10^ F10 ), (0.16243838109804687, F5^ F5  F11^ F11 ), (0.15886083485519145, F5^ F5  F12^ F12 ), (0.15886083485519145, F5^ F5  F13^ F13 ), (0.767389123035518, F6^ F6  F0^ F0 ), (0.767389123035518, F6^ F6  F1^ F1 ), (0.30622126982800785, F6^ F6  F2^ F2 ), (0.30622126982800785, F6^ F6  F3^ F3 ), (0.30208361911523307, F6^ F6  F4^ F4 ), (0.30208361911523307, F6^ F6  F5^ F5 ), (0.2877895677180508, F6^ F6  F6^ F6 ), (0.2877895677180508, F6^ F6  F7^ F7 ), (0.2804819687896131, F6^ F6  F8^ F8 ), (0.2804819687896131, F6^ F6  F9^ F9 ), (0.15475209021245984, F6^ F6  F10^ F10 ), (0.15475209021245984, F6^ F6  F11^ F11 ), (0.15134382699799553, F6^ F6  F12^ F12 ), (0.15134382699799553, F6^ F6  F13^ F13 ), (0.767389123035518, F7^ F7  F0^ F0 ), (0.767389123035518, F7^ F7  F1^ F1 ), (0.30622126982800785, F7^ F7  F2^ F2 ), (0.30622126982800785, F7^ F7  F3^ F3 ), (0.30208361911523307, F7^ F7  F4^ F4 ), (0.30208361911523307, F7^ F7  F5^ F5 ), (0.2877895677180508, F7^ F7  F6^ F6 ), (0.2877895677180508, F7^ F7  F7^ F7 ), (0.2804819687896131, F7^ F7  F8^ F8 ), (0.2804819687896131, F7^ F7  F9^ F9 ), (0.15475209021245984, F7^ F7  F10^ F10 ), (0.15475209021245984, F7^ F7  F11^ F11 ), (0.15134382699799553, F7^ F7  F12^ F12 ), (0.15134382699799553, F7^ F7  F13^ F13 ), (0.7479034551648777, F8^ F8  F0^ F0 ), (0.7479034551648777, F8^ F8  F1^ F1 ), (0.2984456501590826, F8^ F8  F2^ F2 ), (0.2984456501590826, F8^ F8  F3^ F3 ), (0.2944130633377986, F8^ F8  F4^ F4 ), (0.2944130633377986, F8^ F8  F5^ F5 ), (0.2804819687896131, F8^ F8  F6^ F6 ), (0.2804819687896131, F8^ F8  F7^ F7 ), (0.2733599255869175, F8^ F8  F8^ F8 ), (0.2733599255869175, F8^ F8  F9^ F9 ), (0.15082260028140723, F8^ F8  F10^ F10 ), (0.15082260028140723, F8^ F8  F11^ F11 ), (0.147500880233922, F8^ F8  F12^ F12 ), (0.147500880233922, F8^ F8  F13^ F13 ), (0.7479034551648777, F9^ F9  F0^ F0 ), (0.7479034551648777, F9^ F9  F1^ F1 ), (0.2984456501590826, F9^ F9  F2^ F2 ), (0.2984456501590826, F9^ F9  F3^ F3 ), (0.2944130633377986, F9^ F9  F4^ F4 ), (0.2944130633377986, F9^ F9  F5^ F5 ), (0.2804819687896131, F9^ F9  F6^ F6 ), (0.2804819687896131, F9^ F9  F7^ F7 ), (0.2733599255869175, F9^ F9  F8^ F8 ), (0.2733599255869175, F9^ F9  F9^ F9 ), (0.15082260028140723, F9^ F9  F10^ F10 ), (0.15082260028140723, F9^ F9  F11^ F11 ), (0.147500880233922, F9^ F9  F12^ F12 ), (0.147500880233922, F9^ F9  F13^ F13 ), (0.41264550253745813, F10^ F10  F0^ F0 ), (0.41264550253745813, F10^ F10  F1^ F1 ), (0.16466330572421778, F10^ F10  F2^ F2 ), (0.16466330572421778, F10^ F10  F3^ F3 ), (0.16243838109804687, F10^ F10  F4^ F4 ), (0.16243838109804687, F10^ F10  F5^ F5 ), (0.15475209021245984, F10^ F10  F6^ F6 ), (0.15475209021245984, F10^ F10  F7^ F7 ), (0.15082260028140723, F10^ F10  F8^ F8 ), (0.15082260028140723, F10^ F10  F9^ F9 ), (0.08321430695009596, F10^ F10  F10^ F10 ), (0.08321430695009596, F10^ F10  F11^ F11 ), (0.08138159334405823, F10^ F10  F12^ F12 ), (0.08138159334405823, F10^ F10  F13^ F13 ), (0.41264550253745813, F11^ F11  F0^ F0 ), (0.41264550253745813, F11^ F11  F1^ F1 ), (0.16466330572421778, F11^ F11  F2^ F2 ), (0.16466330572421778, F11^ F11  F3^ F3 ), (0.16243838109804687, F11^ F11  F4^ F4 ), (0.16243838109804687, F11^ F11  F5^ F5 ), (0.15475209021245984, F11^ F11  F6^ F6 ), (0.15475209021245984, F11^ F11  F7^ F7 ), (0.15082260028140723, F11^ F11  F8^ F8 ), (0.15082260028140723, F11^ F11  F9^ F9 ), (0.08321430695009596, F11^ F11  F10^ F10 ), (0.08321430695009596, F11^ F11  F11^ F11 ), (0.08138159334405823, F11^ F11  F12^ F12 ), (0.08138159334405823, F11^ F11  F13^ F13 ), (0.4035573895111223, F12^ F12  F0^ F0 ), (0.4035573895111223, F12^ F12  F1^ F1 ), (0.16103675769568102, F12^ F12  F2^ F2 ), (0.16103675769568102, F12^ F12  F3^ F3 ), (0.15886083485519145, F12^ F12  F4^ F4 ), (0.15886083485519145, F12^ F12  F5^ F5 ), (0.15134382699799553, F12^ F12  F6^ F6 ), (0.15134382699799553, F12^ F12  F7^ F7 ), (0.147500880233922, F12^ F12  F8^ F8 ), (0.147500880233922, F12^ F12  F9^ F9 ), (0.08138159334405823, F12^ F12  F10^ F10 ), (0.08138159334405823, F12^ F12  F11^ F11 ), (0.07958924346013586, F12^ F12  F12^ F12 ), (0.07958924346013586, F12^ F12  F13^ F13 ), (0.4035573895111223, F13^ F13  F0^ F0 ), (0.4035573895111223, F13^ F13  F1^ F1 ), (0.16103675769568102, F13^ F13  F2^ F2 ), (0.16103675769568102, F13^ F13  F3^ F3 ), (0.15886083485519145, F13^ F13  F4^ F4 ), (0.15886083485519145, F13^ F13  F5^ F5 ), (0.15134382699799553, F13^ F13  F6^ F6 ), (0.15134382699799553, F13^ F13  F7^ F7 ), (0.147500880233922, F13^ F13  F8^ F8 ), (0.147500880233922, F13^ F13  F9^ F9 ), (0.08138159334405823, F13^ F13  F10^ F10 ), (0.08138159334405823, F13^ F13  F11^ F11 ), (0.07958924346013586, F13^ F13  F12^ F12 ), (0.07958924346013586, F13^ F13  F13^ F13 )
First two body rotation:  [[-0.991 -0.09  -0.    -0.    -0.094  0.    -0.01 ]
 [ 0.109 -0.641 -0.    -0.    -0.587 -0.     0.483]
 [-0.     0.    -0.694  0.     0.    -0.72  -0.   ]
 [ 0.044  0.399  0.    -0.    -0.798  0.    -0.449]
 [ 0.    -0.    -0.    -1.     0.     0.    -0.   ]
 [ 0.057 -0.65  -0.     0.     0.101  0.    -0.751]
 [-0.     0.    -0.72   0.     0.     0.694  0.   ]]

Note

The corresponding rotation matrix for the constant energy term is the identity.

Circuit representations of the rotation operators may be computed using the inquanto.ansatzes.restricted_basis_rotation_to_circuit() method, which uses the RealRestrictedBasisRotationAnsatz to encode the rotations in the Jordan-Wigner picture (see here for more information). For example, below we construct a computable for the expectation value of a single term in (87) with a UCCSD ansatz:

from inquanto.mappings import QubitMappingJordanWigner
from inquanto.ansatzes import FermionSpaceAnsatzUCCSD, CircuitAnsatz, ComposedAnsatz, restricted_basis_rotation_to_circuit
from inquanto.states import FermionState
from inquanto.computables import ExpectationValue

jw = QubitMappingJordanWigner()

operator, rotation = terms[4]  # Selecting the 4th term as an example

ansatz = FermionSpaceAnsatzUCCSD( # |psi>
    fermion_space=7*2,
    fermion_state=FermionState([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]),
    qubit_mapping=jw
)
ansatz_rotator = CircuitAnsatz(restricted_basis_rotation_to_circuit(rotation.T))  # R(U)^
rotated_ansatz = ComposedAnsatz(ansatz_rotator, ansatz)  # R(U)^ |psi>
term_computable = ExpectationValue(rotated_ansatz, jw.operator_map(operator))  # <psi| R(U) X R(U)^ |psi>

Note the transpose of the rotation matrix: \(\hat{R}(\textbf{U}_t)^\dagger = \hat{R}(\textbf{U}_t^T)\), which is true for a real-valued rotation matrix.