I am using the following code in sympy:
from sympy.physics.secondquant import F, Fd, NO, Commutator
from sympy import symbols
a, b, c, d = symbols("a,b,c,d")
comm = NO(Commutator( Fd(a) * F(b), Fd(c) * F(d) ).doit())
print(comm)
# gives me 0
Clearly, this should not be zero.
Well, maybe I do not understand sympy. Here is what I want to calculate:
[F†_a F_b, F†_c F_d] with not necessarily equal indices a, b, c and d.