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$, are bounded uniformly over $a_1,\dots,a_n$.

By permuting we may assume that $R_1 \leq \dots \leq R_n$.  By rescaling we can also normalize the set $\{a_1,\dots,a_n\}$ to have diameter $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).$$

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 $x$, 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.  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 \leq 1$ and $R_n \sim 1$ (and there are only $O(\log^{n-2} (1/R_1))$ choices of the remaining scales $R_2,\dots,R_{n-1}$ for any fixed $R_1,R_n$).