The series of squares of harmonic numbers
$$
\sum_{n=1}^\infty H_n^2 \tag1
$$ is divergent since $\displaystyle \lim_{n \to \infty} H_n^2 \ne 0$, actually from the classic result (6.3.18),
$$
H_n=\ln n+\gamma+\frac1{2n}+O\left(\frac1{n^2}\right),\qquad \, n \to \infty,
$$ where $\gamma$ is the Euler-Mascheroni constant, one gets as $n \to \infty$,
$$
H_n^2=\left(\ln n+\gamma+\frac1{2n} \right)^2+O\left(\frac{\ln n}{n^2}\right)\tag2
$$which tends to infinity.
Then the following new series
$$
\sum_{n=1}^\infty \color{grey}{\left[\color{#151515}{\: H_n^2-\left(\ln n+\gamma+\frac1{2n} \right)^2}\: \right]} \tag3
$$
may be seen as a sort of regularization of $(1)$. The series $(3)$ is absolutely convergent as one may directly deduce from the comparison test with a Bertrand series, using $(2)$.
Question. What is a closed form of $(3)$?
My intuition is that $(3)$ admits a closed form in terms of known constants (or here). I've used the advanced Inverse Symbolic Calculator ISC $2.0$ but it found nothing. My recent attempt, not yet fruitful, has been to convert $(3)$ into an integral representation, starting with
$$
\begin{align}
-\int_0^1\!\left(\!\frac1{\ln x}+\frac1{1-x}-\frac12\!\right)\!x^{n-1}\:dx&=H_n-\ln n-\gamma-\frac1{2n},\quad n\ge1,
\end{align}
$$ and trying to employ similar techniques used here.
Analogous considerations are here or here, one may also explore some variations of $(3)$, like
$$
\sum_{n=1}^\infty \color{grey}{\left[\color{#151515}{\: H_n^q-\left(\ln n+\gamma+\frac1{2n} \right)^q}\: \right]}, \,\sum_{n=1}^\infty (-1)^n \!\color{grey}{\left[\color{#151515}{\: H_n^q-\left(\ln n+\gamma+\frac1{2n} \right)^q}\: \right]}.
$$
Best Answer
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.