The last exercise can be solved using the symbolic Python package called SymPy. SymPy is a Python package for general purpose symbolic algebra. There is a physics module with several interesting submodules. Among these, the submodule called secondquant, contains several functionalities that allow us to test our algebraic manipulations using Wick's theorem and operators for second quantization.

from sympy import *
from sympy.physics.secondquant import *

i, j = symbols('i,j', below_fermi=True)
a, b = symbols('a,b', above_fermi=True)
p, q = symbols('p,q')
print simplify(wicks(Fd(i)*F(a)*Fd(p)*F(q)*Fd(b)*F(j), keep_only_fully_contracted=True))

The code defines single-particle states above and below the Fermi level, in addition to the genereal symbols \( pq \) which can refer to any type of state below or above the Fermi level. Wick's theorem is implemented between the creation and annihilation operators Fd and F, respectively. Using the simplify option, one can lump together several Kronecker-\( \delta \) functions.