Probability that a random pair of points are opposite corners of a square in an $n\times n$ integer lattice
Solution 1:
For given $k$, $1\leq k\leq n-1$, an axis-aligned square of side length $n-k$ can be placed in $k^2$ ways, and each such square hosts $n-k$ squares with vertices on its rim. The total number $Q_n$ of admissible squares therefore is given by $$Q_n=\sum_{k=1}^{n-1}k^2(n-k)=\ldots={n^2(n^2-1)\over12}\ .$$ I think the exact form of the final result is a coincidence, cf. the formula for $\sum_{k=1}^n k^3$. – Working with the individual choices is made more difficult by the following fact: The chosen points have to fulfill a parity condition before you even can think of a possible square. If $n$ is odd the probability that this parity condition holds is $\ne{1\over2}$.
By the way: Say "a random pair of lattice points" since two points chosen independently can coincide.
Solution 2:
I think what you outline looks like a fairly optimal (and elegant) solution. As you pointed out, once we have the count of squares inside the grid, the claim is an immediate consequence, by counting in two ways the number of triples $(x,y,Q)$, where $[x,y]$ is a pair of points and $Q$ is a square, and $xy$ is a diagonal of $Q$.
The formula $6(\# Q) ={n^2\choose 2}$ also has a simple proof: We count separately the number of squares that have $(p,q)$ as their bottom edge, with $p\ge 1$, $q\ge 0$. Since this fits into a square of side length $p+q$ with sides parallel to the axes (and this is the smallest such square), we have $(n-p-q)^2$ such squares. This gives $$ \# Q = \sum_{p\ge 1, q\ge 0} (n-p-q)^2 = \sum_{j=1}^n\sum_{p+q=j} (n-j)^2 = \sum_{j=1}^n j(n-j)^2 = \sum_{k=0}^{n-1} k^2(n-k) , $$ as claimed.
While this answer probably just reproduces what you already did yourself, I don't think it's realistic too expect an alternative solution that's much simpler still.