I guess the answer to this question depends on how we define the rate of convergence. However if we focus on the limiting behavior of characteristic function then the answer is pretty easy. Define by $Z_n := \frac{S_n-2 n}{\sqrt{n \log(n)}}$ as the detrended and normalized random variable and denote by $\kappa_{Z_n}(k)$ the characteristic function of that variable. Note that the characteristic function of $Z_n$ always exists however the higher moments -- as we will see- don't. We have:
\begin{eqnarray}
\kappa_{Z_n}(k) &:=& \int\limits_{[0,1]^n} e^{\imath k \left(
\frac{\sum\limits_{j=1}^n \frac{1}{\sqrt{x_j}} - 2 n}{\sqrt{n \log(n)}}\right)} \prod\limits_{j=1}^n d x_j\\
&=& e^{-\imath \frac{k 2 n}{\sqrt{n \log(n)}}} \cdot \left( e^{\imath \frac{k }{\sqrt{n \log(n)}}}(1+\frac{\imath k}{\sqrt{n \log(n)}}) + \frac{k^2}{n \log(n)}(-\imath \pi +Ei(\frac{\imath k}{\sqrt{n \log(n)}})\right)^n
\end{eqnarray}

Now, we have:
\begin{eqnarray}
\kappa_{Z_n}(k) - e^{-k^2/2} = \frac{1+2 \gamma-\imath \pi +2 \log(k) - \log(\log(n))}{2 \log(n)} \cdot k^2 + O(k^3 \log(k))
\end{eqnarray}

Therefore for "small" values of $k$ we can see from the above what the rate of convergence is. If we want to be more precise and speak about convergence in terms of some norm more work is required. Here I only plot the difference between the characteristic function in question and the limiting function as a function of $k$. The coloring of the curves corresponds is like in the rainbow with violet and red representing $n=2$ and $n=10$ respectively.

**Update:** It is only now that I got familiar with Berry-Esseen-type of bounds by looking at Terence Tao's blog post https://terrytao.wordpress.com/2010/01/05/254a-notes-2-the-central-limit-theorem/#more-3281 . Let us therefore try to make my answer more rigorous. As stated in the question we are looking to find estimates on the supremum norm of the cumulative distribution functions. Let us denote $G := N(0,1)$ and $\phi(x):= 1_{x>0}$ and write:
\begin{align}
\lVert \Psi^{<}_{Z_n} - \Psi^{<}_G \rVert &=
\mbox{sup}_{a\in {\mathbb R}} \left| P(Z_n < a) - P(G < a)\right| \\
&=\mbox{sup}_{a\in {\mathbb R}} \left| E[ \phi(Z_n-a) ] - E[\phi(G-a)]\right| \\
&=\mbox{sup}_{a\in {\mathbb R}} \frac{1}{2\pi} \left|\int\limits_{{\mathbb R}} \hat{\phi}(k) \cdot e^{\imath k a}\cdot \left( \kappa_{Z_n}(k) - e^{-\frac{k^2}{2}}\right) dk\right|\\
&\leq\frac{1}{2\pi} \left|\int\limits_{{\mathbb R}} \left|\hat{\phi}(k)\right| \cdot 1 \cdot \left| \kappa_{Z_n}(k) - e^{-\frac{k^2}{2}}\right| dk\right|
\end{align}
where $\hat{\phi}(k) := \int\limits_{\mathbb R} \phi(x) \exp(\imath k x) dx$. As explained in Terence's Tao blog post $\left| \hat{\phi}(k) \right|$ is bounded for large values of $k$.

Now, all we need to do is to find some upper bound on the log-characteristic function of $Z_n$ for large values of $n$. We have:
\begin{align}
\log(\kappa_{Z_n}(k))
&= -\imath \frac{k 2 n}{\sqrt{n \log(n)}}+n \log\left(
e^{\imath \frac{k }{\sqrt{n \log(n)}}}(1+\frac{\imath k}{\sqrt{n \log(n)}}) + \frac{k^2}{n \log(n)}(-\imath \pi +Ei(\frac{\imath k}{\sqrt{n \log(n)}}) \right)\\
&=-\frac{k^2}{2} + \frac{k^2}{2} \cdot \frac{1+2 \gamma-\imath \pi +2 \log(k) - \log(\log(n))}{\log(n)} \\
&\qquad+ \frac{(\imath k)^3}{3} \cdot \frac{-2+6 \gamma -3 \imath \pi+6 \log(k) - 3 \log(n)-3 \log(\log(n))}{\sqrt{n} [\log(n)]^{3/2}} + O(\frac{k^4 [\log(k)]^2}{n})
\end{align}
where in the last line I used Mathematica's `Series[]`

command to expand the log in a series to the fourth order.
Now clearly we have:
\begin{align}
\left| \kappa_{Z_n}(k) - e^{-\frac{k^2}{2}}\right| &=
e^{-\frac{k^2}{2}} \cdot \left| e^{\frac{k^2}{2} \cdot \frac{1+2 \gamma-\imath \pi +2 \log(k) - \log(\log(n))}{\log(n)}+\frac{(\imath k)^3}{3} \cdot \frac{-2+6 \gamma -3 \imath \pi+6 \log(k) - 3 \log(n)-3 \log(\log(n))}{\sqrt{n} [\log(n)]^{3/2}} + O(\frac{k^4 [\log(k)]^2}{n})} - 1\right|\\
&\leq e^{-\frac{k^2}{2}} \cdot \left(\frac{k^2}{2} \cdot |\frac{1+2 \gamma-\imath \pi +2 \log(k) - \log(\log(n))}{\log(n)}|+ O\!\left(\frac{k^3 [\log(k)]^1}{\sqrt{n \log(n)}} \right)\right)
\end{align}
To summarize we have the following:
\begin{equation}
|| \Psi^{<}_{Z_n} - \Psi^{<}_G ||\le \frac{1}{2\pi} \cdot \int\limits_{{\mathbb R}} |\hat{\phi}(k)| \cdot e^{-\frac{k^2}{2}} \cdot \frac{k^2}{2} dk \cdot \frac{\log(\log(n))}{\log(n)} + O\left(\frac{1}{\log(n)}\right)
\end{equation}
as $n$ goes to infinity.

**Update':** As we can see from the above the rate of convergence tends to be very slow. Of course it may be that our upper bound is not tight and there exists a more accurate upper bound that diminishes much faster with $n$.
For the time being I am not quite sure how to check whether this is the case or not. However some insight can be achieved by doing a Monte Carlo simulation. Below I plot the sample cdf (Blue) along with the limit cdf (Purple) for different values of $n$. Here I took $n=10^1,\cdots,10^4$ and in each case I used $m=10000$ instances of random variables. I used the following Mathematica code to run the simulation:

```
SetOptions[ListPlot, ImageSize -> 500,
LabelStyle -> {15, FontFamily -> "Arial"},
BaseStyle -> {15, FontFamily -> "Bold"}];
SetOptions[ListLogLogPlot, ImageSize -> 500,
LabelStyle -> {15, FontFamily -> "Arial"},
BaseStyle -> {15, FontFamily -> "Bold"}];
ns = Floor[10^Array[# &, 100, {1, 4}]];
m = 10000; myList = {};
Do[
n = ns[[i]];
X = RandomReal[{0, 1}, {m, n}];
ll = (Total[1/Sqrt[#]] & /@ X - 2 n)/Sqrt[n Log[n]];
delta = 1/10;
bins = Table[-5 + delta/2 + j delta, {j, 1, (10 - delta)/delta}];
emp = EmpiricalDistribution[ll];
DD = CDF[emp, bins];
limD = CDF[NormalDistribution[0, 1], bins];
myList =
Join[myList, {n, Max[#],
First[bins[[Ordering[#, 1, #1 > #2 &]]]]} & /@ {Abs[
DD - limD]}];
pl = ListPlot[Transpose[{bins, #}] & /@ {DD, limD}, Joined :> True,
PlotMarkers -> Automatic, PlotLabel -> "n=" <> ToString[n]];
If[MemberQ[{10, 100, 1000, 10000}, n],
Export["CDFs_n_" <> ToString[n] <> ".jpg", pl, "JPEG"]];
PrintTemporary["n=", n, "done"];
, {i, 1, Length[ns]}];
MatrixForm[myList];
pl1 = GraphicsGrid[{{ListLogLogPlot[{#[[1]], #[[2]]} & /@ myList],
ListPlot[{#[[1]], #[[3]]} & /@ myList]}}];
Export["LimitBehavior.jpg", pl1, "JPEG"];
Import["LimitBehavior.jpg"]
```

Here are the plots:

CDFs at n=10

CDFS at n=100

CDFs at n=1000

CDFs at n=10000

Below we also show both the Kolmogorow distance and the location of the maximum as a function of the sample size (left and right respectively).

Kolmogorow distance and location of maximum

As above for $m=500000$ realizations.

Conclusions:

- My Monte Carlo simulation matches the results given by Clement C, despite the smaller number of realizations which I took to avoid running out of memory.
Interestingly enough the location of the maximum seems to be corresponding to $x=0$ and to be independent of the sample size(it would be interesting to find some theoretical explanation of that fact).

With this number of samples($m=10000$) it is hard to tell whether the line in the double logarithmic plot is a straight line or instead concave or convex and more simulations are required to give a definite answer. However if we increase the number of samples to $m=500000$-- which we can achieve by averaging the relevant piece of code over an addition fifty realizations -- we clearly end up with a line that is convex (see the bottom figure above). This rules out a power-law behavior.