This won't help at all in terms of the theory, but a possible approximation approach is to discretize so that the operators are all matrices, and then get your approximated operator by solving a linear algebra problem.
The related problem ``what operators can be represented as a commutator" has being extensively studied for bounded operators. Say, in finite dimensional space the solution exists if and only if the trace of the right-hand side operator is zero
K. Shoda, Einige Satze uber Matrizen, Japan J. Math. 17 (1936), 361-365.
In the paper
Brown, Arlen; Pearcy, Carl. Structure of commutators of operators. Ann. Math. (2) 82, 112-127 (1965)
you can find a complete description of operators in a Hilbert space that are representable as commutator of two bounded linear operators. There are many papers that deal with this problem in other Banach spaces.
In some of these papers it is mentioned that the same problem for unbounded operators appears somewhere in quantum mechanics, but with this I cannot help you because I know nothing about the subject.
It seems to me that a crucial point about your problem is that you don't just want to realize a given operator as a commutator, but as the commutator with a fixed (given) operator. While I don't know the answer to your question, you could try to attack it by studying the action on eigenfunctions for D. If [D,A]=iaLz, then for any eigenfunction f of D, you get a differential equation on Af . Maybe this helps in solving the problem or spotting problems.
I see another possible approach, provided that one has a good understanding of [D,Lz] and of iterated commutators of this form. If this is the case, one could try to understand the Lie algebra generated by D and Lz .
If unknown operator A is a pseudo-differential operator you can try to solve your equation by using the formula for the commutator of pseudo-differential operators.
Basically I'm trying to find a unitary operator U=exp(-iA) such that
U^dagger P^2/2m U = P^2/2m - Omega Lz,
U^dagger Lz U = Lz - mr^2 Omega,
U^dagger V(r,phi,z) U = V(r,phi+Omega t,z),
where V(r,phi,z) is a potential.
The operator U=exp(-i Omega/hbar (t*Lz+mr^2*phi) does the job on Lz and V(r,phi,z), but it doesn't transform the Laplacian as I expect. There are some residual terms in the transformed Laplacian.