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.