How to use the FermionOperator class

The FermionOperator class is used to represent a linear combination of products of fermionic creation and annihilation operators. For example:

\[0.5 a_{0, \alpha}^\dagger a_{3, \alpha} - 0.25 a_{3, \alpha}^\dagger a_{0, \alpha} + (1 + i) a_{1, \beta}^\dagger a_{5, \beta} a_{4, \alpha}^\dagger.\]

Such operators are represented programmatically as a dictionary (hash table) that maps each term to its associated coefficient. An individual term is represented as a tuple of fermionic actions. Each action is described by three pieces of information:

1. Whether the action is to create or destroy a fermion.

2. The spin of the orbital being acted upon (alpha or beta).

3. The numerical index of the orbital.

Internally, the action is itself represented as a tuple containing these three pieces of information, but we recommend using the helper functions cre_a, cre_b, des_a, and des_b to construct these tuples.

The following code example shows how to construct the operator shown above.

import ffsim

op1 = ffsim.FermionOperator(
    {
        (ffsim.cre_a(0), ffsim.des_a(3)): 0.5,
        (ffsim.cre_a(3), ffsim.des_a(0)): -0.25,
        (ffsim.cre_b(1), ffsim.des_b(5), ffsim.cre_a(4)): 1 + 1j,
    }
)
op1

FermionOperator({
    (cre_a(0), des_a(3)): 0.5,
    (cre_a(3), des_a(0)): -0.25,
    (cre_b(1), des_b(5), cre_a(4)): 1+1j
})

Use repr to view a string representation that displays the action tuples.

repr(op1)

'FermionOperator({((True, False, 0), (False, False, 3)): 0.5+0j, ((True, False, 3), (False, False, 0)): -0.25+0j, ((True, True, 1), (False, True, 5), (True, False, 4)): 1+1j})'

FermionOperators support arithmetic operations. Note that when multiplying a FermionOperator by a scalar, the scalar must go on the left, i.e. 2 * op and not op * 2.

op2 = ffsim.FermionOperator(
    {
        (ffsim.cre_b(2),): 1j,
        (ffsim.des_a(3), ffsim.des_b(3)): -0.25,
    }
)

op3 = op1 * op2 + 2 * op1 - op2 / 4

op3

FermionOperator({
    (cre_a(3), des_a(0), cre_b(2)): 0-0.25j,
    (des_a(3), des_b(3)): 0.0625,
    (cre_a(3), des_a(0)): -0.5,
    (cre_b(1), des_b(5), cre_a(4), des_a(3), des_b(3)): -0.25-0.25j,
    (cre_b(1), des_b(5), cre_a(4), cre_b(2)): -1+1j,
    (cre_a(0), des_a(3)): 1,
    (cre_b(1), des_b(5), cre_a(4)): 2+2j,
    (cre_a(0), des_a(3), des_a(3), des_b(3)): -0.125,
    (cre_b(2)): 0-0.25j,
    (cre_a(0), des_a(3), cre_b(2)): 0+0.5j,
    (cre_a(3), des_a(0), des_a(3), des_b(3)): 0.0625
})

It is good to be aware that some in-place operations are especially efficient because they avoid copying data into a new FermionOperator object. These operations are:

• In-place addition and subtraction

• In-place multiplication and division by a scalar

Some examples of usage:

op3 += op1
op3 -= op2
op3 *= 4
op3 /= 1j
op3

FermionOperator({
    (cre_a(3), des_a(0), cre_b(2)): -1,
    (des_a(3), des_b(3)): 0-1.25j,
    (cre_a(3), des_a(0)): 0+3j,
    (cre_b(1), des_b(5), cre_a(4), des_a(3), des_b(3)): -1+1j,
    (cre_b(1), des_b(5), cre_a(4), cre_b(2)): 4+4j,
    (cre_a(0), des_a(3)): 0-6j,
    (cre_b(1), des_b(5), cre_a(4)): 12-12j,
    (cre_a(0), des_a(3), des_a(3), des_b(3)): 0+0.5j,
    (cre_b(2)): -5,
    (cre_a(0), des_a(3), cre_b(2)): 2,
    (cre_a(3), des_a(0), des_a(3), des_b(3)): 0-0.25j
})

Operators can be normal-ordered by calling the normal_ordered method, which returns a new FermionOperator. In the normal-ordered form of a FermionOperator, the operators comprising each term appear from left to right in descending lexicographic order by (action, spin, orb). That is, all creation operators appear before all annihilation operators; within creation/annihilation operators, spin beta operators appear before spin alpha operators, and larger orbital indices appear before smaller orbital indices.

op3.normal_ordered()

FermionOperator({
    (cre_a(0), des_a(3)): 0-6j,
    (cre_b(2), cre_b(1), cre_a(4), des_b(5)): 4+4j,
    (cre_a(3), des_b(3), des_a(3), des_a(0)): 0+0.25j,
    (cre_b(1), cre_a(4), des_b(5), des_b(3), des_a(3)): -1+1j,
    (cre_b(2), cre_a(3), des_a(0)): -1,
    (cre_b(2)): -5,
    (cre_a(3), des_a(0)): 0+3j,
    (des_b(3), des_a(3)): 0+1.25j,
    (cre_b(2), cre_a(0), des_a(3)): 2,
    (cre_b(1), cre_a(4), des_b(5)): -12+12j
})

If a FermionOperator conserves particle number and the Z component of spin, then it can be converted to a LinearOperator using the ffsim.linear_operator function. In order for a FermionOperator to conserve particle number and the Z component of spin, it must preserve the number of spin-alpha fermions and the number of spin-beta fermions. So far, the operators we have created do not satisfy these criteria:

print(op3.conserves_particle_number())
print(op3.conserves_spin_z())

False
False

Let’s construct a FermionOperator that does conserve particle number and spin, and then convert it to a LinearOperator within the subspace of 3 spatial orbitals with 1 alpha and 1 beta fermion.

op4 = ffsim.FermionOperator(
    {
        (ffsim.cre_a(0), ffsim.des_a(3)): 1 + 1j,
        (ffsim.cre_a(1), ffsim.des_a(1), ffsim.cre_b(1), ffsim.des_b(1)): -0.5,
    }
)
print(op4.conserves_particle_number())
print(op4.conserves_spin_z())

norb = 3
nelec = (1, 1)
linop = ffsim.linear_operator(op4, norb=norb, nelec=nelec)

True
True

A LinearOperator can be matrix-multiplied onto a vector:

vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec))
linop @ vec

array([0.        +0.j        , 0.        +0.j        ,
       0.        +0.j        , 0.        +0.j        ,
       0.00380274+0.03040046j, 0.        +0.j        ,
       0.        +0.j        , 0.        +0.j        ,
       0.        +0.j        ])

It can also be passed into most linear algebra routines in scipy.sparse.linalg.

import scipy.sparse.linalg

eigs, vecs = scipy.sparse.linalg.eigs(linop, which="LM", k=1)
eigs

array([-0.5+1.21984743e-18j])