Let's consider this simple example: a particle in a one-dimensional parabolic potential. In the following, $p$ is the particle momentum operator, $m$ is the mass, $x$ is the position operator, and $\omega$ is related to the potential depth.
$$ \begin{aligned} \left[x, p \right] &= xp - px = i \hbar \\ H &= \frac{p^{2}}{2 m} + \frac{1}{2}m \omega^{2} x^{2} \\ a &= \frac{1}{\sqrt{2 m \hbar \omega}} (m \omega x + i p) \end{aligned} $$ $H$ is the Hamiltonian and $a$ is the creation operator.
Let say I would like to evaluate $\left[ H, a \right]$.
from sympy import *
from sympy.physics.quantum import Commutator, Dagger, Operator
m, omega, h = symbols("m, \omega, \hbar", real=True)
x = Operator('x')
p = Operator('p')
H = (p**2 / m + m * omega**2 * x**2) / 2
a = 1 / sqrt(2 * m * h * omega) * (m * omega * x + I * p)
comm = Commutator(H, a)
simplify(expand(comm.doit()))
This snippet of code will produce:
$$ \frac{\sqrt{2} \omega \left(- i \omega m p x^{2} + i \omega m x^{2} p + p^{2} x - x p^{2}\right)}{4 \sqrt{\hbar \omega m}} $$
At this point, I would like to make sympy aware of the fact that $\left[x, p \right] = xp - px = i \hbar$, and then simplifying it further to achieve the final result ($\left[ H, a \right] = -\frac{\hbar \omega}{2} a$). How can I do it?