Calculating commutators in quantum mechanics symbolically with the help of Mathematica

I think I've figured it out, but I haven't checked it super carefully so buyer beware:

Unprotect[NonCommutativeMultiply];
ClearAll[NonCommutativeMultiply];
NonCommutativeMultiply[] := 1;
NonCommutativeMultiply[a_] := a;
NonCommutativeMultiply[first___, const_?NumericQ*b_, rest___] := 
  const*NonCommutativeMultiply[first, b, rest];
NonCommutativeMultiply[first___, const_?NumericQ, rest___] := 
  const*NonCommutativeMultiply[first, rest];

MakeBoxes[
   NonCommutativeMultiply[first___, 
    args : Longest@Repeated[x_, {2, \[Infinity]}], rest___], form_] :=

  
  RowBox[Flatten@{
     If[Length[{first}] > 0,
      {MakeBoxes[NonCommutativeMultiply[first], form], "**"},
      Nothing
      ],
     SuperscriptBox[MakeBoxes[x, form], ToBoxes[Length[{args}], form]],
     If[Length[{rest}] > 0,
      {"**", MakeBoxes[NonCommutativeMultiply[rest], form]},
      Nothing
      ]
     }];
MakeBoxes[NonCommutativeMultiply[arg_], form_] := MakeBoxes[arg, form]

SetAttributes[NonCommutativeMultiply, Flat];
q[i_] ** p[i_] := I + p[i] ** q[i];
q[i_] ** p[j_] := p[j] ** q[i];
q[i_] ** q[j_] /; ! OrderedQ[{i, j}] := q[j] ** q[i]
p[i_] ** p[j_] /; ! OrderedQ[{i, j}] := p[j] ** p[i]

a_ ** (b_ + c_) := a ** b + a ** c;
(b_ + c_) ** a_ := b ** a + c ** a;

(* Allowing for powers in input and output *)
p /: p[i_]^n_Integer := 
  NonCommutativeMultiply @@ ConstantArray[p[i], n];
q /: q[i_]^n_Integer := 
  NonCommutativeMultiply @@ ConstantArray[q[i], n];

As you can see, I went with p and q to denote the conjugate variables to distinguish them from xyz as the directions. You will need to use ** to multiply them together, but you can use powers as well:

comm[a_, b_] := a ** b - b ** a
comm[q[x], p[x]]
comm[q[x], p[y]]
comm[p[x] ** q[x], p[x]]
comm[q[x] ** p[x]^3, p[x]]

I - 2 p[x] ** q[x]

0

-2 p[x]^2 ** q[x] + I p[x]

-I p[x]^3 + 2 p[x]^4 ** q[x]

Edit

You can add the following definition if you want to be able to raise to symbolic powers:

q[i_] ** p[i_]^n_ := I n p[i]^(n - 1) + p[i]^n ** q[i]

James F. Feagin's Quantum Methods with Mathematica book has an elegant implementation of this in chapter 15.1 Commutator Algebra.

It's along the lines of @Sjoerd's answer (but figured I'd provide the reference to the book above), first defining typical identities for the NonCommutativeMultiply symbol:

Unprotect[NonCommutativeMultiply];

A_ ** (B_ + C_) := A ** B + A ** C
(B_ + C_) ** A_ := B ** A + C ** A
A_ ** c_?NumberQ := c A
c_?NumberQ ** A_ := c A
A_ ** (B_ c_?NumberQ) := c A ** B
(A_ c_?NumberQ) ** B_ := c A ** B
A_ ** (B_ c_Rational) := c A ** B
(A_ c_Rational) ** B_ := c A ** B
A_ ** (B_ c_Power) := c A ** B
(A_ c_Power) ** B_ := c A ** B

and then defining a fundamental commutation expression. e.g. for OP's case:

commutator[A_, B_] := A ** B - B ** A
fundamentalCommutation[expr_] := ExpandAll[expr //. p[i_] ** q[i_] :> q[i] ** p[i] - I h]

which indeed recovers the derivative action of the momentum operator:

h /: NumberQ[h] = True;
{commutator[p[x]/(-I h), q[x]], 
commutator[p[x]/(-I h), q[x] ** q[x]],
commutator[p[x]/(-I h), q[x] ** q[x] ** q[x]]} //fundamentalCommutation

{1, 2 q[x], 3 q[x] ** q[x]}

It is then easy to use a different fundamental commutation expression, e.g. for working with raising and lowering operators:

fundamentalComm[expr_] := ExpandAll[expr //. a ** ad :> ad ** a + 1]
{commutator[a, ad], commutator[ad, a], commutator[a ** ad, ad]} // fundamentalComm

{1, -1, ad}

That can be done with the NCAlgebra (non-commutative algebra) package, see the documentation.

Example:

(* Import package *)
<< NC`
<< NCAlgebra`
<< NCGBX`

SetNonCommutative[x, y, px, py]
SetMonomialOrder[x, y, px, py] (* x to the left, p to the right *)
NCSetOutput[NonCommutativeMultiply -> True] (* pretty output *)

(* commutation relations *)
gb = NCMakeGB[{
   x ** y - y ** x,
   px ** py - py ** px,
   x ** py - py ** x,
   y ** px - px ** y,
   x ** px - px ** x - I,
   y ** py - py ** y - I}, 20];

(* Define expression that should be simplified *)
Comm[A_, B_] := A ** B - B ** A;
H = px ** px + py ** py + x ** x + y ** y + x ** x ** x ** x + y ** y ** y ** y;
expression = Comm[H, px] // NCExpand;

NCReplaceRepeated[expression, gb] // NCExpand
--> I x - x ** px ** x + x ** x ** px + I x ** x ** x - x ** px ** x ** x ** x + x ** x ** x ** x ** px

According to the docs, I think the above should have worked, but apparently it doesn't. It seems NCReplaceRepeated and NCExpand have to be applied multiple times until the result "converges":

foo[X_] := NCExpand@NCReplaceRepeated[X, gb]
expression // foo // foo // foo // foo
--> 2 I x + 4 I x ** x ** x