How to sample from a distribution with submodular weights
There is a paper by Gotovos, Hassani, and Krause (NIPS, 2015) titled "Sampling from Probabilistic Submodular Models" which describes a MCMC method based on Gibbs sampling. They show that the mixing time $t_{\text{mix}}(\epsilon)$ is bounded above by 
$$ t_{\text{mix}}(\epsilon) \leq 2 n^2 \exp(\zeta_f) \log \left( \frac{1}{\epsilon p_\min} \right) $$
where $p_\min = \min_{S} p(S)$ and the quantity $\zeta_f$ is
$$\zeta_f =  \max_{A, B \subset V} \left| f(A) + f(B) - f(A \cup B) - f(A \cap B) \right| $$
which I don't know of appearing elsewhere in the literature. There is another paper by Rebeschini and Karbasi (JMLR, 2015) titled "Fast Mixing for Discrete Point Processes" which proposes a different Gibbs sampler whose mixing time is bounded using the so-called "discrete hessian" which I also haven't seen used elsewhere in the literature.

It seems that submodularity of $f$ alone isn't enough to get fast mixing times for MCMC methods. Current methods certainly need to consider other quantities of $f$ to ensure fast mixing times.