Question: How can we evaluate $$\sum_{n=1}^\infty\frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n},$$where $H_n=\frac11+\frac12+\cdots+\frac1n$?
Quick Results
This series converges because $$\frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n}=O\left(\frac{\ln^2n}{n^{3/2}}\right).$$
My Attempt
Recall the integral representation of harmonic number $$H_n=\int_0^1\frac{1-x^n}{1-x}d x$$
we have $$
S=\sum_{n=1}^\infty\frac1n\frac{\binom{2n}n}{4^n}\iint_{[0,1]^2}\frac{(1-x^n)(1-y^n)}{(1-x)(1-y)}d xd y\\
=\tiny\iint_{[0,1]^2}\frac{x y \log (4)-2 x y \log \left(\sqrt{1-x}+1\right)-2 x y \log \left(\sqrt{1-y}+1\right)+2 x y \log \left(\frac{1}{2} \left(\sqrt{1-x y}+1\right)\right)}{\left(\sqrt{1-x y}-1\right) \left(\sqrt{1-x y}+1\right)}dxdy\\
$$
This integral is too hard for me and Mathematica to compute. Numerical integration returns $12.6178$, it agrees with the numerical summation of the original series. I tried to integrate with respect to $x$, but failed.
This is not a complete solution but some first steps.
EDIT 12.04.19 23:20
Much simpler single integral derived.
Original post
The sum in question is
$$s = \sum_{n=1}^\infty a_n\tag{1}$$
with
$$a_n = \frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n}\tag{2}$$
1. Representation as a single integral
1.1
Let us replace just one harmonic number in $a_n$.
Using the definition
$$H_n = \sum _{k=1}^{\infty } \frac{n}{k (k+n)}\tag{3}$$
and writing
$$\frac{1}{n+k}=\int_0^1 x^{n+k-1}\,dx\tag{4}$$
gives for the n-sum
$$\sum_{n=1}^{\infty } \frac{\binom{2 n}{n} H_n x^n}{4^n}=\frac{\partial}{\partial{c}} \left( {_2}F{_1} \left( \frac{1}{2},1,c,x\right)\right)|_{ c \to 1}\tag{5}$$
The remaining k-sum is easily done
$$-\sum _{k=1}^{\infty } \frac{x^{k-1}}{k} =\frac{\log (1-x)}{x} $$
Hence $s$ can be expressed as
$$s_1 = \frac{\partial}{\partial{c}} i(c)|_{ c \to 1} \tag{6a}$$
with
$$i(c) = \int_0^1 \frac{\log (1-x)}{x} {_2}F{_1} \left( \frac{1}{2},1,c,x\right)\,dx\tag{6b}$$
Here ${_2}F{_1}$ is the hypergeometric function.
Numerically, we find in this form
$$s = 12.6216...$$.
1.2 Simpler single integral
The expression derived in the previous paragraph is correct but not very useful because it contains the hypergeometric function. Here we derive the following simpler formula with an elementary integrand.
$$s_2 = \int_0^\infty \frac{v}{\sinh \left(\frac{v}{2}\right)} \left(\frac{v}{\sqrt{2-e^{-v}}}-2 \log \left(\frac{\sqrt{2-e^{-v}}+1}{e^{-\frac{v}{2}}+1}\right)\right)\,dv\tag{7}$$
This is a well converging integral, suited for numerical evaluation. The integrand is depicted here
The derivation starts with replacing both $H_n$ by (3) and (4).
This gives the integral
$$s = \int_0^1 \int_0^1 \frac{\log(1-x) \log(1-y)}{2(1-x y )^{\frac{3}{2}}}\,dx\,dy\tag{8}$$
Transforming $x\to 1-e^{-u}$, $y\to 1-e^{-v}$ leads to
$$s = \int_0^\infty \int_0^v (u v ) \frac{e^{\frac{u+v}{2}}}{(e^u + e^v -1 )^{\frac{3}{2}}}\,du\,dv\tag{8}$$
Here we have made use of the symmetry of the integrand to restrict the integration region to $u\le v$ (and applying a factor 2). Luckily the $u$-integral can be done with the result (7).
2. Sum with asymptotic summands
An attempt to get a feeling for the ingredients of a possible closed form.
The leading asymptotic term of $a_n$ is
$$a_n \simeq b_n = \frac{(\log (n)+\gamma )^2}{\sqrt{\pi } n^{3/2}}\tag{1}$$
The sum of $b_n$ instead of $a_n$ gives
$$s \simeq \sum_{n=1}^\infty b_n = \frac{1}{\sqrt{\pi }}\left(\zeta ''\left(\frac{3}{2}\right)-2 \gamma \zeta '\left(\frac{3}{2}\right)+\gamma ^2 \zeta \left(\frac{3}{2}\right)\right)\simeq 12.0733\tag{2}$$
Here $\zeta(x)$ is the Riemann zeta function (and it derivatives), and $\gamma$ is the Euler-Mascheroni constant.
Notice that the numerical value is close to the one mentioned in the OP. Taking higher terms in the asymptotic expansion of $a_n$ leads to slightly higher numerical values.