Quantifying how crowded points are in $d$ dimensional cube Let $x_1, \cdots, x_n$ be distinct points in the $d$ dimensional cube $[-1,1]^d$. What is a lower bound on the quantity:
$$ \sum_{1 \le j < k \le n} \frac{1}{||x_k-x_j||}$$
where $|| \cdot ||$ is the Euclidean norm.
An asymptotic answer for large $n$ is also fine and bounds on any quantity that looks similar, such as replacing the norm with norm squared, is also welcomed. 
My motivation is a problem that appeared in the book The Cauchy-Schwarz Master Class which proved the following:

If $-1 \le x_1 < x_2 < \cdots < x_n \le 1$ then $$  \sum_{1 \le j < k \le n} \frac{1}{x_k-x_j}  \ge \frac{n^2 \log n}8.$$


(Update below)
The situation seems to be less interesting asymptotically for fixed $d \geq 2$. We have 
$$\frac{1}{||x-y||_2} \geq \frac{1}{2\sqrt{d}}$$ for any $x, y \in [-1, 1]^d$, so there's an easy lower bound of 
$$\sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||} \geq \frac{1}{2\sqrt{d}} \binom{n}{2} = \Omega(n^2).$$ This lower bound is actually achieved (up to a constant) by letting the $x_j$ lie on a grid, say $\Lambda^+ = \{n^{-1/d}, 2n^{-1/d}, \dots, 1\}^d$. 
First note that if we set 
$$\Lambda = \{-1, \dots, -2n^{-1/d}, -n^{-1/d}, 0, n^{-1/d}, 2n^{-1/d}, \dots, 1\}^d$$
we have 
$$\sum_{j \neq k} \frac{1}{||x_k - x_j||} < \sum_{v \in \Lambda, v \neq 0} \frac{1}{||v||}$$
for each $k$. 
Also, since $d \geq 2$, the integral 
$$\int_{[-1, 1]^d} \frac{dx}{||x||}$$
actually converges, say to some $I_d < \infty$, meaning in particular that 
$$\frac{1}{n} \sum_{v \in \Lambda, v \neq 0} \frac{1}{||v||} \to I_d$$
(this is just a standard grid approximation).
Putting these together, we get
$$\sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||} = \frac{1}{2} \sum_k \sum_{j \neq k} \frac{1}{||x_k - x_j||} < \frac{n}{2} \sum_{v \in \Lambda, v \neq 0} \frac{1}{||v||} \sim \frac{n^2}{2} I_d$$
so our sum is $O(n^2)$, matching the lower bound.
Update:
It's possible to prove more. Fix some positive integer $d$ and real $p > 0$. Then we can show that for $x_1, \dots, x_n \in [-1, 1]^d$, the following tight asymptotic bounds hold: when $p < d$, we have
$$ \sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||^p} \geq \Omega(n^2)$$
and when $p > d$,
$$ \sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||^p} \geq \Omega(n^{1 + \frac{p}{d}})$$
and both of these bounds are achieved by taking $\{x_1, \dots, x_n\} = \Lambda^+$. So far I haven't been able to find a tight lower bound for the case $p = d$.
The proof for the case $p < d$ is the same as the proof above, using the fact that the integral 
$$\int_{[-1, 1]^d} \frac{dx}{||x||^p}$$ 
converges.
To prove the lower bound for $p > d$: surround each $x_k$ with a cube $C_k = x_k + [-5n^{-1/d}, 5n^{-1/d}]^d$ of side length $10n^{-1/d}$. If at least $n/2$ of these cubes $C_k$ were disjoint from all other $C_j$, then these $n/2$ cubes, all disjoint, would have a total volume of $10^d/2$, but would all be contained in $[-1-5n^{-1/d}, 1+5n^{-1/d}]$, which has volume at most $4^d$ for sufficiently large $n$, a contradiction. This means we have at least $n/2$ cubes, each of which intersects some other cube, so we get at least $n/2$ ordered pairs $(x_k, x_j)$, each with $||x_k - x_j|| < 10 n^{-1/d} \sqrt{d}$, hence at least $n/4$ such pairs with $j < k$, so
$$\sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||^p} \geq \frac{n}{4}\left(\frac{n^{p/d}}{10^p d^{p/2}}\right) = \Omega(n^{1 + p/d}).$$
To evaluate at $\{x_1, \dots, x_n\} = \Lambda^+$, note that for any $k$, 
$$\sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||^p} < n^{p/d} \cdot \sum_{v \in \mathbb{Z}^d, v \neq 0} \frac{1}{||v||^p}$$
but the sum on the RHS converges to some $S_{d, p} < \infty$, since the integral 
$$\int_{x \in \mathbb{R}^d, ||x|| \geq 1} \frac{dx}{||x||^p}$$
converges, so
$$\sum_{1 \leq j < k \leq n} \frac{1}{||x_k - x_j||^p} = \frac{1}{2} \sum_k \sum_{j \neq k} < \frac{n}{2} \cdot n^{p/d} S_{d, p} = O(n^{1 + p/d}).$$
In the case $p = d$, evaluating the sum on $\Lambda^+$ gives a result of $\Theta(n^2 \log n)$, but I don't know how to show this is optimal for $d > 1$.
