By dyadic decomposition, it suffices to obtain bounds on the quantities
$$
\int_{|x-a_i| \sim R_i \forall i} \left[ \frac{V(x,a_1,\dots,a_n)}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}}\ dx$$
which, when summed over dyadic powers of two $R_1,\dots,R_n > 0$, is bounded by $O(1)$, uniformly over $a_1,\dots,a_n$.

By permuting we may assume that $R_1 \leq \dots \leq R_n$.  The most important geometric feature of the set $\{a_1,\dots,a_n\}$ turns out to be its diameter; by rescaling we can normalize this diameter to be $1$.  The triangle inequality then lets us assume $R_n \gtrsim 1$.  In fact the triangle inequality leaves us with just two scenarios: the "far field" scenario
$$ R_1 \sim \dots \sim R_n \gtrsim 1 \quad (1)$$
and the "near field" scenario
$$ R_1 \leq \dots \leq R_n \sim 1 \quad (2).$$

One could refine the range of possibilities of the near field further by inspecting the finer metric geometry of the $a_1,\dots,a_n$ beyond just using the diameter normalization, but fortunately we do not have to do so for this problem.

As implicitly observed by the OP, the volume $V(x,a_1,\dots,a_n)$ is of size at most $O(R_1 \cdots R_n)$ by multiplying all the lengths emenating from $x$.  On the other hand, since $a_1,\dots,a_n$ lie in a diameter one set and $x$ is at distance $O(R_1)$ from $a_1$, we also have the bound $O(R_1)$ by the base-times-height formula. [Note that this is a significantly superior bound in the far field case, which was identified by the OP as the case where the previous estimates were poor.]  Finally, the condition $|x-a_1| \sim R_1$ restricts $x$ to a set of measure $O(R_1^n)$.  Thus the above integral can be bounded by
$$ \lesssim \left[\frac{\min(R_1 \cdots R_n, R_1)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n.$$

In the far field case (1), this simplifies to $O(R_n^{-2})$ which is summable to $O(1)$.  In the near field case (2), we bound the above by
$$ \left[\frac{(R_1 \cdots R_n)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n = \prod_{j=1}^n \frac{R_1}{R_j} \leq \frac{R_1}{R_n} \sim R_1$$
which is summable to $O(1)$ since $R_1 \lesssim 1$ (and there are only $O(\log^{n-2} (1/R_1))$ choices of the remaining scales $R_2,\dots,R_n$ for any fixed $R_1$).