There are (n-1)^2 1x1 squares, (n-2)^2 2x2 squares, ..., 1^2=1 nxn square.

Then there are diagonal squares of various slope. Categorizing the slopes by the Delta x and Delta y of the top left sloping side, Delta x can be anywhere from 1 to n-1, with Delta y being anywhere from 1 to n - 1 - DeltaX. The resulting tilted square then occupies DeltaX + DeltaY rows and columns, as the upper right sloping side of the square interchanges deltaX and deltaY, and can therefore be treated like an orthogonal square of side DeltaX + DeltaY.

So we get Sigma{r=1 to n-1}(r^2) + Sigma{dx=1 to n-1}Sigma{dy = 1 to n - 1 - dx}((n - (dx+dy))^2)

n         # of squares

1             0
2             1
3             6
4             20
5             50
6             105
7             196
8             336
9             540
10            825

DEFDBL A-Z
FOR n = 1 TO 10
 tot = 0
 FOR r = 1 TO n - 1
   tot = tot + r * r
 NEXT
 FOR dx = 1 TO n - 1
   FOR dy = 1 TO n - 1 - dx
       r = dx + dy
       tot = tot + (n - r) * (n - r)
   NEXT dy
 NEXT dx
 PRINT n, tot
NEXT

The On-Line Encyclopedia of Integer Sequences (Sloane) lists this as sequence A002415. One of its many definitions is

 a(n) = number of squares with corners on an n X n grid. See also A024206, A108279.
 
and its closed-form formula is

n^2*(n^2-1)/12

and they are called the 4-dimensional pyramidal numbers.

Actually, the program works just as well using:

DEFDBL A-Z
FOR n = 1 TO 10
 tot = 0
 FOR dx = 0 TO n - 1
   FOR dy = 1 TO n - 1 - dx
       r = dx + dy
       tot = tot + (n - r) * (n - r)
   NEXT dy
 NEXT dx
 PRINT n, tot
NEXT

corresponding to

Sigma{dx=0 to n-1}Sigma{dy = 1 to n - 1 - dx}((n - (dx+dy))^2)