The similar integral
$$
\int_0^1\frac{\log^2(1-x)}x\mathrm dx=2\zeta(3)
$$
is evaluated in this blog post using the substitution $u=-\log(1-x)$:
$$
\begin{align}
\int_0^1\frac{\log^2(1-x)}x\mathrm dx
&=
\int_0^\infty\frac{u^2}{1-\mathrm e^{-u}}\mathrm e^{-u}\,\mathrm du
\\
&=
\int_0^\infty u^2\sum_{n=1}^\infty\mathrm e^{-nu}\mathrm du
\\
&=
\sum_{n=1}^\infty\int_0^\infty u^2\mathrm e^{-nu}\mathrm du
\\
&=
\sum_{n=1}^\infty\frac2{n^3}
\\
&=
2\zeta(3)\;.
\end{align}
$$
Analogously substituting $u=\log(1+x)$ in the present integral leads to an integral up to $\log2$ that can be evaluated in terms of polylogarithms evaluated at $\frac12$:
$$
\begin{align}
&\int_0^{\log2}\frac{\log^2(1+x)}x\mathrm dx
\\
=&
\int_0^{\log2}\frac{u^2}{\mathrm e^u-1}\mathrm e^u\,\mathrm du
\\
=&
\int_0^{\log2}\frac{u^2}{1-\mathrm e^{-u}}\mathrm du
\\
=&
\int_0^{\log2} u^2\sum_{n=0}^\infty\mathrm e^{-nu}\mathrm du
\\
=&
\sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du
\\
=&
\sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du
\\
=&
\frac13\log^32+\sum_{n=1}^\infty\frac1n\left(-2^{-n}\log^22+2\int_0^{\log2} u\mathrm e^{-nu}\mathrm du\right)
\\
=&
\frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22+\frac2{n^2}\left(-2^{-n}\log2+\int_0^{\log2}\mathrm e^{-nu}\mathrm du\right)\right)
\\
=&
\frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22-\frac2{n^2}2^{-n}\log2-\frac2{n^3}\left(2^{-n}-1\right)\right)
\\
=&
\def\Li{\operatorname{Li}}
\frac13\log^32-\Li_1\left(\frac12\right)\log^22-2\Li_2\left(\frac12\right)\log2-2\Li_3\left(\frac12\right)+2\zeta(3)
\\
=&
\frac13\log^32-\log2\log^22-2\left(\frac{\pi^2}{12}-\frac12\log^22\right)\log2-2\left(\frac16\log^32-\frac{\pi^2}{12}\log2+\frac78\zeta(3)\right)+2\zeta(3)
\\
=&
\frac{\zeta(3)}4\;.
\end{align}
$$
Not only is this a rather complicated derivation of a much simpler result; it also looks as if the polylogarithm values may have been obtained using the present integral in the first place.
This is a sequel to my comments to the question which was too long to fit in another comment.
We have the formulas for $\vartheta_{3}^{10}(q), \vartheta_{3}^{12}(q)$ from Topics in Analytic Number Theory by Rademacher (famous for proving an infinite series formula to calculate the number of partitions of a positive integer) on page 198:
\begin{align}
\vartheta_{3}^{10}(q) &= 1 + \frac{4}{5}\left\{\sum_{n = 1}^{\infty}\frac{2n^{4}q^{n}}{1 + q^{2n}} + \sum_{n = 1}^{\infty}(-1)^{n - 1}\frac{(2n - 1)^{4}q^{2n - 1}}{1 - q^{2n - 1}}\right\} + \frac{2}{5}\vartheta_{3}^{2}(q)\vartheta_{2}^{4}(q)\vartheta_{4}^{4}(q)\tag{1}\\
\vartheta_{3}^{12}(q) &= 1 + 8\sum_{n = 1}^{\infty}\frac{n^{5}q^{n}}{1 - q^{2n}} - 8\sum_{n = 1}^{\infty}(-1)^{n}\frac{n^{5}q^{2n}}{1 - q^{2n}} + \vartheta_{2}^{4}(q)\vartheta_{3}^{4}(q)\vartheta_{4}^{4}(q)\tag{2}
\end{align}
Finding a general formula for $\vartheta_{3}^{k}(q)$ for even positive integer $k$ is a difficult problem but using the methods given in Rademacher's book it looks like it is possible to obtain such formulas at the cost of heavy symbolic manipulation for a specific $k$.
Update: I found one pattern in your formulas by using the substitution $x = K'(k)/K(k)$ so that when $x = 0$ then $k = 1$ and when $x = \infty$ then $k = 0$ and moreover $$\frac{dx}{dk} = -\frac{\pi}{2kk'^{2}K^{2}}$$ so that the integral of $\vartheta_{4}^{2n}(e^{-\pi x})/(1 + x^{2})$ is transformed into $$\int_{0}^{1}\left(\frac{2k'K}{\pi}\right)^{n}\frac{1}{K^{2} + K'^{2}}\frac{\pi}{2kk'^{2}}\,dk = \left(\frac{2}{\pi}\right)^{n - 1}\int_{0}^{1}\frac{k^{-1}k^{'(n - 2)}K^{n}}{K^{2} + K'^{2}}\,dk$$ and that explains (at least to some extent) the occurrence of $\dfrac{1}{\pi^{n - 1}}$ in your formulas.
Next it is easy to prove one of the formulas in $(14)$. We have $$\vartheta_{2}^{2}\vartheta_{4}^{2} = kk'(2K/\pi)^{2}$$ and hence $$\int_{0}^{\infty}\vartheta_{2}^{2}\vartheta_{4}^{2}\,dx = \int_{0}^{1}kk'\cdot\frac{4K^{2}}{\pi^{2}}\cdot\frac{\pi}{2kk'^{2}K^{2}}\,dk = \frac{2}{\pi}\int_{0}^{1}\frac{dk}{\sqrt{1 - k^{2}}} = 1$$ I wonder if similar technique can be applied to prove other formulas.
If $q = e^{-\pi x}$ then $dx = -\dfrac{dq}{\pi q}$ and interval $(0, \infty)$ changes to $(0, 1)$ and hence we can express the first integral of $(14)$ as $$\frac{1}{\pi}\int_{0}^{1}\vartheta_{2}^{4}(q)\vartheta_{4}^{2}(q)\,\frac{dq}{q} = \frac{16}{\pi}\int_{0}^{1}\psi^{4}(q^{2})\phi^{2}(-q)\,dq$$ Next $$\psi^{4}(q^{2}) = \sum_{n = 0}^{\infty}\frac{(2n + 1)q^{2n}}{1 - q^{4n + 2}}, \phi^{2}(-q) = 1 + 4\sum_{n = 1}^{\infty}\frac{(-1)^{n}q^{n}}{1 + q^{2n}}$$ I wonder if you can utilize the above Lambert series to prove that the desired integral is equal to $1$. It appears that if we express the integrand as a Lambert series then it can also be expressed as the logarithmic derivative of some product of theta functions and the integral can be evaluated. See this paper regarding some integrals related to theta functions (all of it was given by Ramanujan in his lost notebook).
Best Answer
A solution can be obtained from the formula for $\int_0^\infty\frac{x^x e^{-ax}}{\Gamma(1+x)}\,dx$ given in the linked post (but not proven there), and the asymptotics of Lambert's $\mathrm{W}_{-1}(z)$ as $z\to-1/e$.
In turn, the formula can be obtained from Hankel's integral (also asked for here) $$\frac{1}{\Gamma(s)}=\frac{1}{2\pi i}\int_{-\infty}^{(0^+)}z^{-s}e^z\,dz,$$ where the path of integration encircles the negative real axis (the branch cut of $z^{-s}$).
This is done below (but I prefer to avoid references to "deep properties" of $\mathrm{W}_{-1}$).
For $a>1$ and $x>0$, we put $s=1+x$ and substitute $z=xw$ (the path may be left as is): $$\frac{x^x e^{-ax}}{\Gamma(1+x)}=\frac{1}{2\pi i}\int_{-\infty}^{(0^+)}\left(\frac{e^{w-a}}{w}\right)^x\frac{dw}{w}.$$
Now we deform the contour to have $|e^{w-a}/w|\leqslant c<1$ on it (say, we locate it to the left of $\Re w=b$, and make it enclose $|w|=r$, where $1<r<b<a$). Then integration over $x$ results in an absolutely convergent double integral, and, exchanging the integrations, we obtain ($\lambda$ is the contour described above) $$f(a):=\int_0^\infty\frac{x^x e^{-ax}}{\Gamma(1+x)}\,dx=\frac{1}{2\pi i}\int_\lambda\frac{dw}{w(a-w+\log w)}.$$
The integrand has a pole at $w=w_a$, the solution of $w_a-\log w_a=a$ (here's how Lambert's $\mathrm{W}$ comes in), and the residue is $1/(1-w_a)$; the integral along $|w|=R$ tends to $0$ as $R\to\infty$, giving finally $$f(a)=1/(w_a-1).$$
So, for $a>0$ we have $$\int_0^\infty\left(\frac{x^x e^{-x}}{\Gamma(1+x)}-\frac{1}{\sqrt{2\pi x}}\right)e^{-ax}\,dx=\frac{1}{b}-\frac{1}{\sqrt{2a}},$$ where $b=b(a)$ is the solution of $b-\log(1+b)=a$. The ability to take $a\to 0^+$ under the integral sign is easy to justify (the resulting integral is absolutely convergent, so DCT is applicable).
Thus, the given integral is equal to the limit $$\lim_{b\to 0^+}\left(\frac1b-\frac{1}{\sqrt{2\big(b-\log(1+b)\big)}}\right),$$ which is evaluated easily, and is indeed equal to $-1/3$.