Even if you can derive the logp of a multivariate Bernoulli that behaves as you want, marginalizing by enumeration will likely be computationally prohibitive (since the number of possible terms grows exponentially with size, whereas in the independent cases it’s linear), unless there are some nice symmetries that can be exploited.