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$).