A virial theorem in the Gaussian Unitary Ensemble

In the Gaussian Unitary Ensemble of $n \times n$ matrices, one can calculate the following expectation value $$\begin{equation}\label{eq:20151201a} \mathbb{E} \left[ \sum_{ i \neq j} \dfrac{1}{ \left( \lambda_i - \lambda_j \right)^2} \right] = \dfrac{1}{2} n (n-1) \end{equation}$$ with $\lambda_1 , \ldots, \lambda_n$ the eigenvalues of the random matrix $H$. I have normalized the GUE such that $\mathbb{E} [ H_{ij} H_{kl} ] = \delta_{il} \delta_{jk}$. In this blog post, I check $\eqref{eq:20151201a}$ with a Monte Carlo simulation in Mathematica.
I generate an $n \times n$ random matrix in the GUE with the following code

RandomMatrixGUE[n_] := Module[{A, B, H},
  A = RandomVariate[NormalDistribution[], {n, n}];
  B = RandomVariate[NormalDistribution[], {n, n}];
  H = 1/2 (A + Transpose[A]) + I / 2 (B - Transpose[B]);
  H]

Given a matrix $H$, I calculate the sum $\sum_{ i \neq j} \frac{1}{ \left( \lambda_i - \lambda_j \right)^2}$ with

f[H_] := Module[{lambda = Eigenvalues[H], i, j, n = Length[H], result = 0},
  For[i = 1, i <= n, i++,
  For[j = 1, j < i, j++, result += 1/ (lambda[[i]] - lambda[[j]])^2]];
  2*result] 

I then use
randomsamples[n_, NSamples_] := Table[A = RandomMatrixGUE[n]; f[A], {k, NSamples}]

and

SeedRandom[1]
For [n = 2, n <= 5, n++,
 NSamples = 20000;
 result = randomsamples[n, NSamples];
 vev = Mean[result];
 vev = NumberForm[vev, {3, 2}];
 std = StandardDeviation[result]/Sqrt[NSamples];
 std = NumberForm[std, {3, 2}];
 exact = 1/2 n (n - 1);
 Print["\n n = ", n, ": exact = ", exact, ", Monte Carlo gives ", 
  vev, " with std ", std];
 difference = Abs[vev - exact];
 If[difference <= 3 std, 
  Print["the difference is smaller than 3 std"], 
  Print["The Monte Carlo result is different from exact result, \
namely ", difference/std, " standard deviations"]];]

The last code produces the following output:

 n = 2: exact = 1, Monte Carlo gives 0.94 with std 0.02

 n = 3: exact = 3, Monte Carlo gives 3.20 with std 0.31

 n = 4: exact = 6, Monte Carlo gives 6.06 with std 0.14

 n = 5: exact = 10, Monte Carlo gives 10.30 with std 0.48

The Monte Carlo results above lie within 3 standard deviations of the exact value for $n = 2,3,4,5$. This is thus a good test on $\eqref{eq:20151201a}$.
Remarks:

• The proof of $\eqref{eq:20151201a}$ can be found in the chapter on Brownian motion in Mehta's book. The proof is based on setting up an Ornstein-Uhlenbeck process and then analyzing the related heat equation. The argument is directly taken from a paper by Dyson. I do not know if there is a more direct way of calculating $\eqref{eq:20151201a}$.
• Illustrations of the Dyson Ornstein-Uhlenbeck process can be found in a previous post.