You can use pm.Potential to add arbitrary logp terms to the graph. We discourage that a bit because you lose the 1-to-1 representation between the random and the logp graphs and can’t guarantee correct results when using prior or posterior predictive sampling.
Anyway the Stan model API can be mimicked by something like:
theta = pm.Flat("theta")
pm.Potential("term1", pm.logp(pm.Normal.dist(2, 4), theta))
pm.Potential("term2", pm.logp(pm.Normal.dist(0, 10), theta))
For the likelihood, if I understand correctly you could specify two observed variables with the same data?