FermionOperator

class FermionOperator

Bases: object

A spin-less fermionic operator.

Definition

This operator is defined by a linear combination of products of fermionic creation and annihilation operators acting on spin-less fermionic modes. That is to say, the individual terms fulfill the following anti-commutation relations: [1]

\[\left\{a^\dagger_\alpha, a^\dagger_\beta\right\} = \left\{a_\alpha, a_\beta\right\} = 0,~~\text{and}~~ \left\{a_\alpha, a^\dagger_\beta\right\} = \delta_{\alpha\beta} \, ,\]

where \(\alpha\) and \(\beta\) do not distinguish the spin species of the fermionic modes they are indexing.

This makes the definition of the entire operator the following:

\[\text{\texttt{FermionOperator}} = \sum_i c_i \bigotimes_j \hat{A_j} \, ,\]

where \(\hat{A_j} \in \{ a_j, a^\dagger_j \}\) and \(c_i\) is the (complex) coefficient making up the linear combination of products. The index \(j\) can take any value between 0 and the number of fermionic modes acted upon by the operator minus 1.

Implementation

This class stores the terms and coefficients in multiple sparse vectors, akin to the compressed sparse row format commonly used for sparse matrices. More concretely, a single operator contains 4 arrays:

coeffs

A vector of complex coefficients consisting of two 64-bit floating point numbers.

actions

A vector of booleans storing the nature of the second-quantization actions.

modes

A vector of 32-bit integers storing the fermionic mode indices acted upon.

boundaries

A vector of integers indicating the boundaries in actions and modes.

Entries in actions indicate creation (annihilation) operators by True (False). Fermionic modes indexed by modes are considered spinless.

This data structure allows for very efficient construction and manipulation of operators. However, it implies that duplicate terms may be contained in an operator at any moment. These must be resolved manually through the use of simplify().

Construction

An operator can be constructed directly by providing the arrays outlined above:

>>> from qiskit_fermions.operators import FermionOperator
>>> coeffs = [1.0, 2.0, -3.0, 4.0j, -0.5j]
>>> actions = [True, False, False, True, True, True, False, False]
>>> modes = [0, 0, 0, 1, 0, 1, 2, 3]
>>> boundaries = [0, 0, 1, 2, 4, 8]
>>> op = FermionOperator(coeffs, actions, modes, boundaries)
>>> print(op)
  1.000000e0 +0.000000e0j * ()
 -3.000000e0 +0.000000e0j * (-_0)
  0.000000e0 +4.000000e0j * (-_0 +_1)
  2.000000e0 +0.000000e0j * (+_0)
 -0.000000e0-5.000000e-1j * (+_0 +_1 -_2 -_3)

For convenience, it is possible to construct an operator from a Python dictionary like so:

>>> from qiskit_fermions.operators import cre, ann
>>> op = FermionOperator.from_dict(
...     {
...         (): 1.0,
...         (cre(0),): 2.0,
...         (ann(0),): -3.0,
...         (ann(0), cre(1)): 4.0j,
...         (cre(0), cre(1), ann(2), ann(3)): -0.5j,
...     }
... )
>>> print(op)
  1.000000e0 +0.000000e0j * ()
 -3.000000e0 +0.000000e0j * (-_0)
  0.000000e0 +4.000000e0j * (-_0 +_1)
  2.000000e0 +0.000000e0j * (+_0)
 -0.000000e0-5.000000e-1j * (+_0 +_1 -_2 -_3)

In this example, we have leveraged cre() and ann() for creating the creation and annihilation operators at the specified modes.

In addition, the following construction and quick helper methods are available:

zero()

Constructs the additive identity operator.

one()

Constructs the multiplicative identity operator.

Iteration

Since the underlying data structure is implemented in Rust and has a non-trivial layout, it cannot be iterated over directly:

>>> list(iter(op))
Traceback (most recent call last):
  ...
TypeError: 'qiskit_fermions.operators.fermion_operator.FermionOperator' object is not iterable

Instead, this class provides custom iterators to fulfill this purpose:

>>> list(sorted(op.iter_terms()))
[([], (1+0j)), ([(False, 0)], (-3+0j)), ([(False, 0), (True, 1)], 4j), ([(True, 0)], (2+0j)), ([(True, 0), (True, 1), (False, 2), (False, 3)], (-0-0.5j))]

See also

iter_terms()

For more relevant implementation details.

The table below lists all available iterators:

iter_terms()

An iterator over the operator's terms.

Arithmetics

The following arithmetic operations are supported:

Addition/Subtraction

>>> op = FermionOperator.one()
>>> (op + op).simplify()
FermionOperator.from_dict({(): 2+0j})
>>> (op - op).simplify()
FermionOperator.from_dict({})
>>> op += op
>>> op.simplify()
FermionOperator.from_dict({(): 2+0j})
>>> op -= op
>>> op.simplify()
FermionOperator.from_dict({})

Scalar Multiplication/Divison

>>> op = FermionOperator.one()
>>> (2 * op).simplify()
FermionOperator.from_dict({(): 2+0j})
>>> (op / 2).simplify()
FermionOperator.from_dict({(): 0.5+0j})
>>> op *= 2
>>> op.simplify()
FermionOperator.from_dict({(): 2+0j})
>>> op /= 2
>>> op.simplify()
FermionOperator.from_dict({(): 1+0j})

Operator Composition

Note

Operator composition corresponds to left-multiplication: c = a & b corresponds to \(C = B A\). In other words, the composition of two operators returns a resulting operator that performs “first a and then b”.

>>> op1 = FermionOperator.from_dict({(): 2.0, (cre(0),): 3.0})
>>> op2 = FermionOperator.from_dict({(): 1.5, (ann(1),): 4.0})
>>> comp = (op1 & op2).simplify()
>>> print(comp)
  3.000000e0 +0.000000e0j * ()
  8.000000e0 +0.000000e0j * (-_1)
  1.200000e1 +0.000000e0j * (-_1 +_0)
  4.500000e0 +0.000000e0j * (+_0)
>>> op2 &= op1
>>> print(op2.simplify())
  3.000000e0 +0.000000e0j * ()
  8.000000e0 +0.000000e0j * (-_1)
  4.500000e0 +0.000000e0j * (+_0)
  1.200000e1 +0.000000e0j * (+_0 -_1)
>>> squared = (op1 ** 2).simplify()
>>> print(squared)
  4.000000e0 +0.000000e0j * ()
  1.200000e1 +0.000000e0j * (+_0)
  9.000000e0 +0.000000e0j * (+_0 +_0)

Other Operations

In addition to the magic methods that correspond to the arithmetic operations outlined above, the following methods are available:

adjoint()

Returns the Hermitian conjugate (or adjoint) of this operator.

ichop([atol])

Removes terms whose coefficient magnitude lies below the provided threshold.

simplify([atol])

Returns an equivalent but simplified operator.

normal_ordered()

Returns an equivalent operator with normal ordered terms.

Properties

Finally, various methods exist to check certain properties of an operator:

is_hermitian([atol])

Returns whether this operator is Hermitian.

many_body_order()

Returns the many-body order of this operator.

conserves_particle_number()

Returns whether this operator is particle-number conserving.

Methods

adjoint()

Returns the Hermitian conjugate (or adjoint) of this operator.

This affects the terms and coefficients as follows:

  • the actions in each term reverse their order and flip between creation and annihilation

  • the coefficients are complex conjugated

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): -1.0j, ((True, 0), (False, 1)): 1.0})
>>> adj = op.adjoint()
>>> print(adj)
 -0.000000e0 +1.000000e0j * ()
  1.000000e0 -0.000000e0j * (+_1 -_0)
conserves_particle_number()

Returns whether this operator is particle-number conserving.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({((True, 0), (False, 1)): 1})
>>> op.conserves_particle_number()
True
>>> op = FermionOperator.from_dict({((True, 0),): 1})
>>> op.conserves_particle_number()
False
Returns:

Whether this operator is particle-number conserving.

equiv(other, atol=1e-08)

Checks this operator for equivalence with another operator.

Equivalence in this context means approximate equality up to the specified absolute tolerance. To be more precise, this method returns True, when all the absolute values of the coefficients in the difference other - self are below the specified threshold atol.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): 1e-7})
>>> zero = FermionOperator.zero()
>>> op.equiv(zero)
False
>>> op.equiv(zero, 1e-6)
True
>>> op.equiv(zero, 1e-9)
False
Parameters:

  • other – the other operator to compare with.

  • atol – the absolute tolerance for the comparison. This value defaults to 1e-8.

from_1body_tril_spin(one_body_b, norb)

Constructs an operator from separate spin-species triangular 1-body integrals.

The resulting operator is defined by

\[\sum_i c^\alpha_{ii} a^\dagger_i a_i + c^\beta_{ii} a^\dagger_{i+n} a_{i+n} + \sum_{i \lt j} c^\alpha_{ij} (a^\dagger_i a_j + a^\dagger_j a_i) + c^\beta_{ij} (a^\dagger_{i+n} a_{j+n} + a^\dagger_{j+n} a_{i+n})\]

where \(c^\alpha\) (\(c^\beta\)) are the integral coefficients stored in one_body_a (one_body_b, resp.), \(i\) and \(j\) are the indices expanded from the triangular index \(ij\) which indexes the arrays, and \(n\) is the number of orbitals, norb.

>>> import numpy as np
>>> from qiskit_fermions.operators import FermionOperator
>>> one_body_a = np.array([1.0, 2.0, 3.0])
>>> one_body_b = np.array([-1.0, -2.0, -3.0])
>>> op = FermionOperator.from_1body_tril_spin(one_body_a, one_body_b, norb=2)
>>> print(op)
  1.000000e0 +0.000000e0j * (+_0 -_0)
  2.000000e0 +0.000000e0j * (+_0 -_1)
  2.000000e0 +0.000000e0j * (+_1 -_0)
  3.000000e0 +0.000000e0j * (+_1 -_1)
 -1.000000e0 +0.000000e0j * (+_2 -_2)
 -2.000000e0 +0.000000e0j * (+_2 -_3)
 -2.000000e0 +0.000000e0j * (+_3 -_2)
 -3.000000e0 +0.000000e0j * (+_3 -_3)
Parameters:

  • one_body_a – a 1-dimensional array of length \(n * (n + 1) / 2\) storing the 1-body electronic integral coefficients of the \(\alpha\)-spin species, as a flattened triangular matrix.

  • one_body_b – a 1-dimensional array of length \(n * (n + 1) / 2\) storing the 1-body electronic integral coefficients of the \(\beta\)-spin species, as a flattened triangular matrix.

  • norb – the number of orbitals, \(n\).

Returns:

The 1-body component of the electronic structure Hamiltonian as defined above.

from_1body_tril_spin_sym(norb)

Constructs an operator from spin-symmetric triangular 1-body integrals.

The resulting operator is defined by

\[\sum_i c^\alpha_{ii} (a^\dagger_i a_i + a^\dagger_{i+n} a_{i+n}) + \sum_{i \lt j} c^\alpha_{ij} (a^\dagger_i a_j + a^\dagger_j a_i + a^\dagger_{i+n} a_{j+n} + a^\dagger_{j+n} a_{i+n})\]

where \(c^\alpha\) are the integral coefficients stored in one_body_a, \(i\) and \(j\) are the indices expanded from the triangular index \(ij\) which indexes the array, and \(n\) is the number of orbitals, norb.

>>> import numpy as np
>>> from qiskit_fermions.operators import FermionOperator
>>> one_body_a = np.array([1.0, 2.0, 3.0])
>>> op = FermionOperator.from_1body_tril_spin_sym(one_body_a, norb=2)
>>> print(op)
  1.000000e0 +0.000000e0j * (+_0 -_0)
  2.000000e0 +0.000000e0j * (+_0 -_1)
  2.000000e0 +0.000000e0j * (+_1 -_0)
  3.000000e0 +0.000000e0j * (+_1 -_1)
  1.000000e0 +0.000000e0j * (+_2 -_2)
  2.000000e0 +0.000000e0j * (+_2 -_3)
  2.000000e0 +0.000000e0j * (+_3 -_2)
  3.000000e0 +0.000000e0j * (+_3 -_3)
Parameters:

  • one_body_a – a 1-dimensional array of length \(n * (n + 1) / 2\) storing the 1-body electronic integral coefficients of the \(\alpha\)-spin species, as a flattened triangular matrix.

  • norb – the number of orbitals, \(n\).

Returns:

The 1-body component of the electronic structure Hamiltonian as defined above.

from_2body_tril_spin(two_body_ab, two_body_bb, norb)

Constructs an operator from separate spin-species triangular 2-body integrals.

The resulting operator is defined by

\[\sum_{ijkl} \frac{1}{2} \sum_{(i,j,k,l) \in \mathcal{P}(ijkl)} c^{\alpha\alpha}_{ijkl} a^\dagger_i a^\dagger_k a_l a_j + c^{\beta\beta}_{ijkl} a^\dagger_{i+n} a^\dagger_{k+n} a_{l+n} a_{j+n} + \sum_{ijkl} \frac{1}{2} \sum_{(i,j,k,l) \in \mathcal{P'}(ijkl)} c^{\alpha\beta}_{ijkl} a^\dagger_{i+n} a^\dagger_k a_l a_{j+n} + c^{\alpha\beta}_{ijkl} a^\dagger_i a^\dagger_{k+n} a_{l+n} a_j +\]

where \(c^{\alpha\alpha}\) (\(c^{\alpha\beta}\), \(c^{\beta\beta}\)) are the integral coefficients stored in two_body_aa (two_body_ab, two_body_bb, resp.), \(ijkl\) is the running index of the array, \(\mathcal{P}\) (\(\mathcal{P'}\)) generates the unique permutations of the 4-index \((i,j,k,l)\) (see below), and \(n\) is the number of orbitals, norb.

Note

two_body_aa and two_body_bb are a S8-fold symmetric arrays. That means, they are the flattened lower-triangular data of matrices of shape (npair, npair), where npair = (norb * (norb + 1) // 2. These in turn are the lower-triangular data of the 4-dimensional arrays of shape (norb, norb, norb, norb). Therefore, \(\mathcal{P}\) above expands the flattened index \(ijkl\) into all index permutations \((i,j,k,l)\) that index these 4-dimensional arrays.

However, two_body_ab is only S4-fold symmetric. Thus, it contains the full data of the (npair, npair) matrix (but still in flattened form). \(\mathcal{P'}\) performs the corresponding index expansion. (In the definition above, we reused the index \(ijkl\) as an abuse of notation.)

>>> import numpy as np
>>> from qiskit_fermions.operators import FermionOperator
>>> two_body_aa = np.arange(1, 7, dtype=float)
>>> two_body_ab = np.arange(11, 20, dtype=float)
>>> two_body_bb = np.arange(-1, -7, -1, dtype=float)
>>> op = FermionOperator.from_2body_tril_spin(two_body_aa, two_body_ab, two_body_bb, norb=2)
>>> len(op)
64
Parameters:

  • two_body_aa – a 1-dimensional array of the S8-fold symmetric 2-body electronic integral coefficients of the \(\alpha\alpha\)-spin species, as a flattened array.

  • two_body_ab – a 1-dimensional array of the S4-fold symmetric 2-body electronic integral coefficients of the \(\alpha\beta\)-spin species, as a flattened array.

  • two_body_bb – a 1-dimensional array of the S8-fold symmetric 2-body electronic integral coefficients of the \(\beta\beta\)-spin species, as a flattened array.

  • norb – the number of orbitals, \(n\).

Returns:

The 2-body component of the electronic structure Hamiltonian as defined above.

from_2body_tril_spin_sym(norb)

Constructs an operator from spin-symmetric triangular 2-body integrals.

The resulting operator is defined by

\[\sum_{ijkl} \frac{1}{2} c^{\alpha\alpha}_{ijkl} \sum_{(i,j,k,l) \in \mathcal{P}(ijkl)} (a^\dagger_i a^\dagger_k a_l a_j + a^\dagger_{i+n} a^\dagger_k a_l a_{j+n} + a^\dagger_i a^\dagger_{k+n} a_{l+n} a_j + a^\dagger_{i+n} a^\dagger_{k+n} a_{l+n} a_{j+n})\]

where \(c^{\alpha\alpha}\) are the integral coefficients stored in two_body_aa, \(ijkl\) is the running index of the array, \(\mathcal{P}\) generates the unique permutations of the 4-index \((i,j,k,l)\) (see below), and \(n\) is the number of orbitals, norb.

Note

two_body_aa is an S8-fold symmetric array. That means, it is the flattened lower-triangular data of a matrix of shape (npair, npair), where npair = (norb * (norb + 1) // 2. This in turn is the lower-triangular data of the 4-dimensional array of shape (norb, norb, norb, norb). Therefore, \(\mathcal{P}\) above expands the flattened index \(ijkl\) into all index permutations \((i,j,k,l)\) that index this 4-dimensional array.

>>> import numpy as np
>>> from qiskit_fermions.operators import FermionOperator
>>> two_body_aa = np.arange(1, 7, dtype=float)
>>> op = FermionOperator.from_2body_tril_spin_sym(two_body_aa, norb=2)
>>> len(op)
64
Parameters:

  • two_body_aa – a 1-dimensional array of the S8-fold symmetric 2-body electronic integral coefficients of the \(\alpha\alpha\)-spin species, as a flattened array.

  • norb – the number of orbitals, \(n\).

Returns:

The 2-body component of the electronic structure Hamiltonian as defined above.

from_dict()

Constructs a new operator from a dictionary.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict(
...     {
...         (): 1.0-1.0j,
...         ((True, 0), (False, 1)): 2.0,
...     }
... )
>>> print(op)
  1.000000e0 -1.000000e0j * ()
  2.000000e0 +0.000000e0j * (+_0 -_1)
Parameters:

data – a dictionary mapping tuples of terms to complex coefficients. Each key is a tuple of (bool, int) pairs. You may use cre() and ann() to simplify their construction.

Returns:

A new operator.

from_fcidump()

Constructs a FermionOperator from an FCIDump data structure.

Assuming you have an FCIDump file called molecule.fcidump, you can construct the second-quantized operator like so:

from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.library import FCIDump

fcidump = FCIDump.from_file("molecule.fcidump")
operator = FermionOperator.from_fcidump(fcidump)
Parameters:

fcidump – the FCIDump data structure.

Returns:

The constructed operator.

ichop(atol=1e-08)

Removes terms whose coefficient magnitude lies below the provided threshold.

Caution

This method truncates coefficients greedily! If the acted upon operator may contain separate coefficients for duplicate terms consider calling simplify() instead!

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): 1e-4, ((True, 0),): 1e-6, ((False, 0),): 1e-10})
>>> print(op)
  1.000000e-4 +0.000000e0j * ()
 1.000000e-10 +0.000000e0j * (-_0)
  1.000000e-6 +0.000000e0j * (+_0)
>>> op.ichop()
>>> print(op)
  1.000000e-4 +0.000000e0j * ()
  1.000000e-6 +0.000000e0j * (+_0)
>>> op.ichop(1e-5)
>>> print(op)
  1.000000e-4 +0.000000e0j * ()
Parameters:

atol – the absolute tolerance for the cutoff. This value defaults to 1e-8.

is_hermitian(atol=1e-08)

Returns whether this operator is Hermitian.

Note

This check is implemented using equiv() on the normal_ordered() difference of self and its adjoint() and zero().

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({
...     ((True, 0), (False, 1)): 1.00001j,
...     ((True, 1), (False, 0)): -1j,
... })
>>> op.is_hermitian()
False
>>> op.is_hermitian(1e-4)
True
Parameters:

atol – The numerical accuracy upto which coefficients are considered equal. This value defaults to 1e-8.

Returns:

Whether this operator is Hermitian.

iter_terms()

An iterator over the operator’s terms.

Warning

Mutating the iteration items does not affect the underlying operator data.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): 2.0, ((True, 0),): 1.0, ((False, 1),): -1.0j})
>>> list(sorted(op.iter_terms()))
[([], (2+0j)), ([(False, 1)], (-0-1j)), ([(True, 0)], (1+0j))]
many_body_order()

Returns the many-body order of this operator.

Note

The many-body order is defined as the length of the longest term contained in the operator.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({
...     ((True, 0), (False, 1), (True, 2), (False, 3)): 1,
... })
>>> op.many_body_order()
4
Returns:

The many-body order of this operator.

normal_ordered()

Returns an equivalent operator with normal ordered terms.

The normal order of an operator term is defined such that all creation actions before all annihilation actions and the modes of actions within each group descend lexicographically (e.g. +_1 +_0 -_1 -_0).

Note

When a term is being reordered, the anti-commutation relations have to be taken into account, \(a_i a^\dagger_j = \delta_{ij} - a^\dagger_j a^i\), implying that the number of terms may change.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({((False, 1), (True, 1), (False, 0), (True, 0)): 1})
>>> print(op.normal_ordered().simplify())
  1.000000e0 +0.000000e0j * ()
 -1.000000e0 +0.000000e0j * (+_0 -_0)
 -1.000000e0 +0.000000e0j * (+_1 -_1)
 -1.000000e0 +0.000000e0j * (+_1 +_0 -_1 -_0)
Returns:

An equivalent but normal-ordered operator.

one()

Constructs the multiplicative identity operator.

Composing the operator that is constructed by this method with another one has no effect.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): 2.0})
>>> one = FermionOperator.one()
>>> op & one == op
True
simplify(atol=1e-08)

Returns an equivalent but simplified operator.

The simplification process first sums all coefficients that belong to equal terms and then only retains those whose total coefficient exceeds the specified tolerance (just like ichop()).

When an operator has been arithmetically manipulated or constructed in a way that does not guarantee unique terms, this method should be called before applying any method that filters numerically small coefficients to avoid loss of information. See the example below which showcases how ichop() can truncate terms that sum to a total coefficient magnitude which should not be truncated:

>>> from qiskit_fermions.operators import FermionOperator
>>> coeffs = [1e-5] * int(1e5)
>>> boundaries = [0] + [0] * int(1e5)
>>> op = FermionOperator(coeffs, [], [], boundaries)
>>> canon = op.simplify(1e-4)
>>> assert canon.equiv(op.one(), 1e-6)
>>> op.ichop(1e-4)
>>> assert op.equiv(op.zero(), 1e-6)
Parameters:

atol – the absolute tolerance for the cutoff. This value defaults to 1e-8.

Returns:

An equivalent but simplified operator.

zero()

Constructs the additive identity operator.

Adding the operator that is constructed by this method to another one has no effect.

>>> from qiskit_fermions.operators import FermionOperator
>>> op = FermionOperator.from_dict({(): 2.0})
>>> zero = FermionOperator.zero()
>>> op + zero == op
True