We can continue along these lines and define a normal-ordered Hamiltonian with respect to a given reference state. In our first step we just define the Hamiltonian
from sympy import symbols, latex, WildFunction, collect, Rational, simplify
from sympy.physics.secondquant import F, Fd, wicks, AntiSymmetricTensor, substitute_dummies, NO, evaluate_deltas
# setup hamiltonian
p,q,r,s = symbols('p q r s',dummy=True)
f = AntiSymmetricTensor('f',(p,),(q,))
pr = Fd(p)*F(q)
v = AntiSymmetricTensor('v',(p,q),(r,s))
pqsr = Fd(p)*Fd(q)*F(s)*F(r)
#define the Hamiltonian
Hamiltonian = f*pr + Rational(1)/Rational(4)*v*pqsr
#define indices for states above and below the Fermi level
index_rule = {
'below': 'kl',
'above': 'cd',
'general': 'pqrs'
}
Hnormal = substitute_dummies(Hamiltonian,new_indices=True, pretty_indices=index_rule)
print "Hamiltonian defined as:", latex(Hnormal)
which results in $$ f^{q}_{p} a^\dagger_{q} a_{p} + \frac{1}{4} v^{sr}_{qp} a^\dagger_{s} a^\dagger_{r} a_{p} a_{q} $$