Here is a nice solution from R. Tauraso, for those interested in this problem:
Solution proposed by Roberto Tauraso, Roma, Italy.
The Multiple Harmonic Sum is defined by $$H_{n}(s_{1}, \dots, s_{l}) == \sum_{0 < k_{1} < k_{2} < \dots < k_{l} \leq n}\prod_{i = 1}^{l}\frac{1}{k_{i}^{s_{i}}}$$ Then, by the known properties of the stuffle product,
\begin{align}
H_{n}^{2}(1) &= 2H_{n}(1, 1) + H_{n}(2)\notag\\
H_{n}^{3}(1) &= 6H_{n}(1, 1, 1) + 3H_{n}(1, 2) + 3H_{n}(2, 1) + H_{n}(3)\notag
\end{align}
Therefore, for $N > 0$,
\begin{align}
\sum_{n = 1}^{N}\frac{H_{n}^{2}}{n} &= 2\sum_{n = 1}^{N}\frac{H_{n}(1, 1)}{n} + \sum_{n = 1}^{N}\frac{H_{n}(2)}{n}\notag\\
&= 2H_{N}(1, 1, 1) + 2H_{N}(1, 2) + H_{N}(3) + H_{N}(2, 1)\notag\\
&= \frac{1}{3}H_{N}^{3}(1) + \frac{2}{3}H_{N}(3) + H_{N}(1, 2)\notag\\
&= \frac{1}{3}H_{N}^{3}(1) + \frac{2}{3}\zeta(3) + \zeta(2, 1) + o(1)\notag\\
&= \frac{1}{3}H_{N}^{3}(1) + \frac{5}{3}\zeta(3) + o(1)\tag{1}
\end{align}
because $\zeta(2, 1) = \zeta(3)$ (see for example Thirty-two Goldbach Variations by J. M. Borwein and D. M. Bradley).
Moreover, $H_{N}(1) = \ln N + \gamma + O(1/N)$ implies that
\begin{align}
\frac{1}{3}H_{N}^{3}(1) - \sum_{n = 1}^{N}\frac{(\gamma + \ln n)^{2}}{n} &= \frac{1}{3}H_{N}^{3}(1) - \gamma^{2}H_{N}(1) - 2\gamma\sum_{n = 1}^{N}\frac{\ln n}{n} - \sum_{n = 1}^{N}\frac{\ln^{2}n}{n}\notag\\
&= \frac{\ln^{3}N}{3} + \frac{\gamma^{3}}{3} + \gamma\ln^{2}N + \gamma^{2}\ln N - \gamma^{2}(\ln N + \gamma)\notag\\
&\,\,\,\,\,\,\,\, -2\gamma\left(\gamma_{1} + \frac{\ln^{2}N}{2}\right) - \left(\gamma_{2} + \frac{\ln^{3}N}{3}\right) + o(1)\notag\\
&= -\frac{2}{3}\gamma^{3} - 2\gamma\gamma_{1} - \gamma_{2} + o(1)\tag{2}
\end{align}
By adding together $(1)$ and $(2)$, and by taking the limit as $N \to \infty$, we get the result.
Just some considerations for now.
I have proved here that
$$ \sum_{n=1}^{N}H_n^2 = (N+1) H_N^2-(2N+1)H_N+2N \tag{$\color{blue}{1}$} $$
and we have:
$$ \mathcal{L}^{-1}\left(\frac{\log x+\gamma}{x}\right)(s)=-\log(s)\tag{2} $$
$$ \mathcal{L}^{-1}\left(\frac{\left(\log x+\gamma\right)^2}{x}\right)(s)=-\zeta(2)+\log^2(s)\tag{3} $$
hence the partial sums of the given series can be rearranged as follows:
$$\begin{eqnarray*}\sum_{n=1}^{N}\left[H_n^2-\left(\log n+\gamma+\frac{1}{2n}\right)^2\right]&=&(\color{blue}{1})-\frac{H_N^{(2)}}{4}-\sum_{n=1}^{N}\frac{\log n+\gamma}{n}-\sum_{n=1}^{N}n\frac{(\log n+\gamma)^2}{n}\\&=&(\color{blue}{1})-\frac{H_N^{(2)}}{4}+\int_{0}^{+\infty}\frac{\log(s)(1-e^{-Ns})}{e^s-1}\,ds\\&-&\int_{0}^{+\infty}\frac{e^{-N s} \left(e^{(1+N) s}+N-e^s (1+N)\right)\left(\zeta(2)-\log^2 s\right)}{\left(-1+e^s\right)^2}\,ds\end{eqnarray*}$$
If we find a way to distribute $(\color{blue}{1})$ over the last two integrals, in a way ensuring they are convergent integrals as $N\to +\infty$, we are done. At first sight the closed form of the LHS appears to be related with $\zeta(2)$, $\zeta'(0)=-\frac{1}{2}\log(2\pi)$ and
$$ \zeta''(0)=\frac{\gamma^2}{2}-\frac{\pi ^2}{24}-\frac{1}{2}\log^2(2\pi)+\gamma_1$$
with $\gamma_1$ being a Stieltjes constant. An alternative, brute-force way is just to compute the asymptotic expansions of
$$ \sum_{n=1}^{N}\log^2(n),\qquad \sum_{n=1}^{N}\frac{\log(n)}{n} $$
with the sufficient degree of accuracy (I guess that to stop at the $O\left(\frac{1}{N^3}\right)$ term is enough), that is just a tedious exercise about summation by parts.
Best Answer
Proof of (*)
Adding the four finite sums,
$$\sum_{k=1}^{n}H_{k}=(n+1)H_{n}-n,$$
$$\sum_{k=1}^{n}\ln{k}=\ln{n!},$$
$$\sum_{k=1}^{n}\gamma=\gamma\,n,$$
$$\sum_{k=1}^{n}\frac{1}{2k}=\frac12H_{n},$$
gives us a representation of the $n$-th partial sum for the infinite series. Writing the infinite series as the limit of partial sums, we get:
$$\begin{align} S &=\sum_{k=1}^{\infty}\left(H_{k}-\ln{n}-\gamma-\frac{1}{2k}\right)\\ &=\lim_{n\to\infty}\sum_{k=1}^{n}\left(H_{k}-\ln{n}-\gamma-\frac{1}{2k}\right)\\ &=\lim_{n\to\infty}\left((n+1)H_{n}-n-\ln{n!}-\gamma\,n-\frac12H_{n}\right)\\ &=\lim_{n\to\infty}\left(\left(n+\frac12\right)H_{n}-(1+\gamma)n-\ln{n!}\right). \end{align}$$
Use Stirling's approximation for the factorial to obtain an asymptotic formula for the log-factorial term in the series:
$$n!\sim\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\\ \implies \ln{n!}\sim\ln{\left(\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\right)}=\left(n+\frac12\right)\ln{n}-n+\frac12\ln{(2\pi)}.$$
Then,
$$\begin{align} S &=\lim_{n\to\infty}\left(\left(n+\frac12\right)H_{n}-(1+\gamma)n-\ln{n!}\right)\\ &=\lim_{n\to\infty}\left(\left(n+\frac12\right)H_{n}-(1+\gamma)n-\left(n+\frac12\right)\ln{n}+n-\frac12\ln{(2\pi)}\right)\\ &=\lim_{n\to\infty}\left(\left(n+\frac12\right)H_{n}-\gamma\,n-\left(n+\frac12\right)\ln{n}\right)-\frac12\ln{(2\pi)}\\ &=\lim_{n\to\infty}\left(n\left(H_{n}-\gamma-\ln{n}\right)+\frac12\left(H_{n}-\ln{n}\right)\right)-\frac12\ln{(2\pi)}\\ &=\lim_{n\to\infty}n\left(H_{n}-\gamma-\ln{n}\right)+\frac12\lim_{n\to\infty}\left(H_{n}-\ln{n}\right)-\frac12\ln{(2\pi)}\\ &=\lim_{n\to\infty}n\left(H_{n}-\gamma-\ln{n}\right)+\frac12\gamma-\frac12\ln{(2\pi)}\\ &=\frac12+\frac12\gamma-\frac12\ln{(2\pi)}.~~~\blacksquare \end{align}$$
Appendix:
Using the asymptotic series for the digamma function given by Eq.16 on this Wolfram Mathworld page,
$$\begin{align} \lim_{n\to\infty}n\left(H_{n}-\gamma-\ln{n}\right) &=\lim_{n\to\infty}n\left(\Psi{(n+1)}-\ln{n}\right)\\ &=\lim_{n\to\infty}n\left(\frac{1}{2n}-\sum_{\ell=1}^{\infty}\frac{B_{2\ell}}{2\ell n^{2\ell}}\right)\\ &=\frac12-\lim_{n\to\infty}\sum_{\ell=1}^{\infty}\frac{B_{2\ell}}{2\ell n^{2\ell-1}}\\ &=\frac12. \end{align}$$