In our next step we define the reference energy \( E_0 \) and redefine the Hamiltonian by subtracting the reference energy and collecting the coefficients for all normal-ordered products (given by the NO function).
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)
E0 = wicks(Hnormal,keep_only_fully_contracted=True)
Hnormal = Hnormal-E0
w = WildFunction('w')
Hnormal = collect(Hnormal, NO(w))
Hnormal = evaluate_deltas(Hnormal)
print latex(Hnormal)
which gives us $$ - f^{i}_{i} + f^{q}_{p} a^\dagger_{q} a_{p} - \frac{1}{4} v^{ii}_{ii} - \frac{1}{4} v^{ii}_{ii} + \frac{1}{4} v^{sr}_{qp} a^\dagger_{r} a^\dagger_{s} a_{q} a_{p}, $$ again as expected, with the reference energy to be subtracted.