A tough integral:$\int_0^{+\infty}\left( \frac1{\log(x+1)-\log x}-x-\frac12\right)^2 dx$ I would like to prove the convergence of

$$I=\int_0^{+\infty}\left( \frac1{\log(x+1)-\log x}-x-\frac12\right)^2 dx$$ 

then obtain a closed form of $I$.
Convergence is ensured by the fact that $x \mapsto  f(x)=\left( \frac1{\log(x+1)-\log x}-x-\frac12\right)^2$ is continuous on $(0,+\infty)$ with $f(x) \sim \dfrac14$ as $x \to 0^+$ and $f(x) \sim \dfrac1{144 x^2}$ as $x \to +\infty.$
I did not succeed in finding a closed form of $I$. My attempt was to consider a certain parameter integral then make some differentiation to get rid of the $\log$ terms...
 A: Here is an alternate solution: We introduce
$$ I(s) = \int_{0}^{\infty} \left( \frac{1}{e^x - 1} + \frac{1}{2} - \frac{1}{x}\right)^2 \frac{x^s e^x}{(e^x - 1)^2} \, dx. $$
It is easy to check $I = I(0)$ and $I(s)$ is analytic for $\Re(s) > -1$. We can also check that
\begin{align*}
J(n, s)
&:= \int_{0}^{\infty} \frac{x^{s-1}e^x}{(e^x - 1)^{n+1}} \, dx
 = \sum_{k=n}^{\infty} \binom{k}{n} \frac{\Gamma(s)}{k^s} \\
&= \frac{\Gamma(s)}{n!} \sum_{j=0}^{n} (-1)^{n-j} \left[{n \atop j}\right] \zeta(s-j), \tag{1}
\end{align*}
where $\left[{n \atop j}\right]$ denotes the Stirling's number of the first kind. Using this, if $\Re(s)$ is large, then we find that
\begin{align*}
I(s) &= J(3, s+1) + J(2,s+1) + \tfrac{1}{4}J(1, s+1) \\
&\qquad - 2J(2, s) - J(1, s) + J(1, s-1). \tag{2}
\end{align*}
Plugging (1) to this equation and using the principle of analytic continuation, we check that the resulting expression is still valid near $s = 0$. Simplifying (2), we get
$$ I(s) = \left( \frac{\Gamma(s+1)}{6} - \Gamma(s) + \Gamma(s-1) \right)\zeta(s-2) + \frac{1}{12}\Gamma(s+1)\zeta(s) $$
and taking $s \to 0$ we have
$$ I(0) = -2\zeta'(-2) + \frac{1}{12}\zeta(0) = \frac{\zeta(3)}{2\pi^2} - \frac{1}{24}. $$
A: We have:
$$ I = \frac{1}{4}\int_{0}^{+\infty}\left(\frac{1}{2}+\frac{1}{e^t-1}-\frac{1}{t}\right)^2\frac{dt}{\sinh^2(t/2)}$$
where:
$$\frac{1}{2}+\frac{1}{e^t-1}-\frac{1}{t} = \sum_{n\geq 1}\frac{2t}{t^2+4\pi^2 n^2}\tag{1}$$
as well as:
$$ \frac{1}{\sinh^2(t)}=\frac{1}{t^2}+2\sum_{n\geq 1}\frac{(t^2-n^2\pi^2)}{(t^2+n^2\pi^2)^2}\tag{2}$$
so our integral depends on the two series:
$$ \sum_{n\geq 1}\sum_{m\geq 1}\int_{0}^{+\infty}\frac{dt}{(t^2+4n^2\pi^2)(t^2+4m^2\pi^2)}=\frac{1}{16\pi^2}\sum_{n\geq 1}\sum_{m\geq 1}\frac{1}{mn(m+n)}=\color{red}{\frac{\zeta(3)}{8\pi^2}}$$
$$ \sum_{c\geq 1}\sum_{a\geq 1}\sum_{b\geq 1}\int_{0}^{+\infty}\frac{t^2(t^2-4c^2\pi^2)\,dt}{(t^2+4\pi^2 a^2)(t^2+4\pi^2 b^2)(t^2+4\pi^2 c^2)^2}\\=\frac{1}{16\pi^2}\sum_{c,a,b\geq 1}\frac{ab-c^2}{(a+b)(a+c)^2(b+c)^2}.\tag{3}$$
The first series is well-known, for instance:
$$\sum_{m,n\geq 1}\frac{1}{mn(m+n)}=\sum_{n\geq 1}\frac{2H_{n-1}}{n^2}=2\int_{0}^{1}\frac{\log(1-x)\log x}{x}\,dx = \int_{0}^{1}\frac{\log^2 x}{1-x^2}\,dx.$$
What about the second one? It can be written as an integral over the unit cube, so I believe it can be computed through the techniques outlined by C. Viola in Birational transformations and values of the Riemann zeta-function, paragraph $4$, Permutation groups for triple integrals. By exploiting the $a\leftrightarrow b$ symmetry we have:
$$\begin{eqnarray*} \sum_{c,a,b\geq 1}\frac{ab-c^2}{(a+b)(a+c)^2(b+c)^2} &=& \sum_{c,a,b\geq 1}\frac{(a+c)(b-c)}{(a+b)(a+c)^2(b+c)^2}\\&=&\sum_{c,a,b\geq 1}\left(\frac{1}{a+c}-\frac{1}{a+b}\right)\frac{1}{(b+c)^2}\end{eqnarray*} $$
suggesting that the last series in $(3)$ is just zero, but there are some convergence issues: we are allowed to exchange $a$ and $b$, but we are not allowed to exchange $b$ and $c$. We have:
$$ \sum_{a,b\geq 1}\left(\frac{1}{a+c}-\frac{1}{a+b}\right)\frac{1}{(b+c)^2}=\frac{1}{c^3}+\frac{\gamma+\psi(c)}{c^2}-\frac{\psi'(c)}{c}-\frac{1}{2}\psi''(c)\tag{4}$$
hence:
$$ \sum_{c\geq 1}\sum_{a,b\geq 1}\frac{ab-c^2}{(a+b)(a+c)^2(b+c)^2}=-\frac{1}{2}\sum_{c\geq 1}\psi''(c)\tag{5}$$
and:

$$ I = \color{red}{\frac{\zeta(3)}{2\pi^2}-\frac{1}{24}} = 0.01923024745\ldots \tag{6}$$

The last step follows from the integral representation for the $\psi$ function:
$$\sum_{c\geq 1}\psi''(c)=-\sum_{s\geq 0}\int_{0}^{1}\frac{x^s\,\log^2 x}{1-x}\,dx = -\int_{0}^{1}\frac{\log^2 x}{(1-x)^2}\,dx = -2\zeta(2).$$
