CustomDist should be fine.
If you are not in a hurry this will actually be trivial after Infer logp for elemwise transformations of multivariate variables by ricardoV94 · Pull Request #6607 · pymc-devs/pymc · GitHub get’s merged, probably in the next release.
In that case all you will need to do is something like
def dist(mu, cov, size):
return pm.math.exp(pm.MvNornal.dist(mu, cov, size=size))
with pm.Model():
mu = ...
cov = ...
x = pm.CustomDist("x", mu, cov, dist=dist)
PyMC will be able to figure out the logp automatically from the dist function. For now you would still implement a logp function.