This is a pure optimization problem. We will center our unit square at the origin, so the particles are located at $(x_i,y_i)$, with $-\tfrac{1}{2}\leq x_i\leq \tfrac{1}{2}$ and $-\tfrac{1}{2}\leq y_i\leq \tfrac{1}{2}$. Our optimization problem is rather simple to state:
\begin{align*}
\underset{x,y}{\text{minimize}}\qquad & \sum_{1\leq i \lt j \leq n} \frac{1}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}} \\
\text{subject to:}\qquad & -\tfrac{1}{2}\leq x_i\leq \tfrac{1}{2} \quad i=1,\dots,n\\
& -\tfrac{1}{2}\leq y_i\leq \tfrac{1}{2} \quad i=1,\dots,n
\end{align*}Taking a numerical approach, I used a local gradient-based solver in Mathematica with randomized initial conditions, solved 50 times, and took the best solution found. I solved the problem for $n=2,\dots,17$ and here were the best solutions found:

Here is a plot of the optimal cost as the number of particles grows:

The case $n=9$
For the case $n=9$, interestingly, the optimal solution is not to place the points in a regular $3\times 3$ grid. Using an irregular arrangement can produce a slightly lower cost, as shown in the figure below.

The optimal arrangement for $9$ particles shown above has the particles at coordinates:
\[
\left(
\begin{array}{cc}
-0.5 & 0.5 \\
0.5 & 0.5 \\
0.5 & -0.5 \\
-0.5 & -0.5 \\
0. & 0.5 \\
0.5 & 0.0469081 \\
-0.5 & 0.0469081 \\
0.181044 & -0.5 \\
-0.181044 & -0.5 \\
\end{array}
\right)
\]
Side note: If we change the energy function to $1/r^2$, then the optimal arrangement does become the regular $3\times 3$ grid.