Preliminaries
Looking at the real part of $\log\left(1-e^{ix}\right)$, we get
$$
\log\left(2\sin\left(\frac x2\right)\right)=-\sum_{k=1}^\infty\frac{\cos(kx)}{k}\tag1
$$
Integrating by parts twice, we get
$$
\int_0^{2\pi}x^2\cos(kx)\,\mathrm{d}x=\left\{\begin{array}{}\frac{8\pi^3}3&\text{if }k=0\\\frac{4\pi}{k^2}&\text{if }k\ne0\end{array}\right.\tag2
$$
In the derivations below, we use the values of $\zeta(2)$ and $\zeta(4)$ computed at the end of this answer.
Start to evaluate one sum, then get another:
$$
\begin{align}
\color{#C00}{\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{(k+j)^2kj}}
&=\sum_{j=1}^\infty\sum_{k=1}^\infty\left(\frac1{kj}-\frac1{k(k+j)}\right)\left(\frac1{kj}-\frac1{j(k+j)}\right)\tag{3a}\\[3pt]
&=\color{#00F}{\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{k^2}\frac1{j^2}}-2\color{#090}{\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{k^2(k+j)j}}+\color{#C00}{\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{(k+j)^2kj}}\tag{3b}\\
\color{#090}{\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{k^2(k+j)j}}
&=\frac12\color{#00F}{\zeta(2)^2}\tag{3c}\\[3pt]
&=\frac{\pi^4}{72}\tag{3d}
\end{align}
$$
Explanation:
$\text{(3a)}$: $\frac1{kj}-\frac1{k(k+j)}=\frac1{j(j+k)}$ and $\frac1{kj}-\frac1{j(k+j)}=\frac1{k(j+k)}$
$\text{(3b)}$: expand the products and separate the sums
$\phantom{\text{(3b):}}$ note that $\sum\limits_{j=1}^\infty\sum\limits_{k=1}^\infty\frac1{k^2(k+j)j}=\sum\limits_{j=1}^\infty\sum\limits_{k=1}^\infty\frac1{k(k+j)j^2}$
$\text{(3c)}$: cancel the red sums, move the green sum to the left side, and divide by $2$
$\text{(3d)}$: simplify
Evaluate the sum we started before (using the other sum we got before):
$$
\begin{align}
\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{(k+j)^2kj}
&=\sum_{j=1}^\infty\sum_{k=j+1}^\infty\frac1{k^2(k-j)j}\tag{4a}\\
&=\sum_{k=1}^\infty\sum_{j=1}^{k-1}\frac1{k^2(k-j)j}\tag{4b}\\
&=\sum_{k=1}^\infty\sum_{j=1}^{k-1}\frac1{k^3}\left(\frac1{k-j}+\frac1j\right)\tag{4c}\\
&=2\sum_{k=1}^\infty\left(\frac{H_k}{k^3}-\frac1{k^4}\right)\tag{4d}\\
&=2\sum_{k=1}^\infty\sum_{j=1}^\infty\frac1{k^3}\left(\frac1j-\frac1{j+k}\right)-2\zeta(4)\tag{4e}\\
&=2\sum_{k=1}^\infty\sum_{j=1}^\infty\frac1{k^2(k+j)j}-2\zeta(4)\tag{4f}\\
&=\frac{\pi^4}{36}-\frac{\pi^4}{45}\tag{4g}\\[6pt]
&=\frac{\pi^4}{180}\tag{4h}
\end{align}
$$
Explanation:
$\text{(4a)}$: substitute $k\mapsto k-j$
$\text{(4b)}$: swap order of summation
$\phantom{\text{(4b):}}$ (we include $k=1$ since the inner sum is then $0$)
$\text{(4c)}$: partial fractions
$\text{(4d)}$: $\sum\limits_{j=1}^{k-1}\frac1{k-j}=\sum\limits_{j=1}^{k-1}\frac1j=H_k-\frac1k$
$\text{(4e)}$: $H_k=\sum\limits_{j=1}^\infty\left(\frac1j-\frac1{j+k}\right)$
$\text{(4f)}$: simplify the summand
$\text{(4g)}$: apply $(3)$
$\text{(4h)}$: simplify
Putting The Preliminaries To Work
$$
\begin{align}
&\int_0^{2\pi}x^2\log\!\left(2\sin\left(\frac{x}2\right)\right)^2\,\mathrm{d}x\\
&=\int_0^{2\pi}x^2\sum_{j=1}^\infty\sum_{k=1}^\infty\frac{\cos(kx)}{k}\frac{\cos(jx)}{j}\,\mathrm{d}x\tag{5a}\\
&=\int_0^{2\pi}x^2\sum_{j=1}^\infty\sum_{k=1}^\infty\frac{\color{#C00}{\cos((k+j)x)}+\color{#090}{\cos((k-j)x)}}{2kj}\,\mathrm{d}x\tag{5b}\\
&=\color{#C00}{2\pi\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{(k+j)^2kj}}+\color{#090}{4\pi\sum_{j=1}^\infty\sum_{k=j+1}^\infty\frac1{(k-j)^2kj}+\frac{4\pi^3}3\sum_{k=1}^\infty\frac1{k^2}}\tag{5c}\\
&=2\pi\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{(k+j)^2kj}+4\pi\sum_{j=1}^\infty\sum_{k=1}^\infty\frac1{k^2(k+j)j}+\frac{4\pi^3}3\zeta(2)\tag{5d}\\
&=\frac{\pi^5}{90}+\frac{\pi^5}{18}+\frac{2\pi^5}9\tag{5e}\\[6pt]
&=\frac{13\pi^5}{45}\tag{5f}
\end{align}
$$
Explanation:
$\text{(5a)}$: apply $(1)$
$\text{(5b)}$: $\cos(a)\cos(b)=\frac{\cos(a+b)+\cos(a-b)}2$
$\text{(5c)}$: apply $(2)$
$\text{(5d)}$: substitute $k\mapsto k+j$ in the middle sum
$\text{(5e)}$: apply $(3)$ and $(4)$
$\text{(5f)}$: simplify
$\newcommand{\+}{^{\dagger}}%
\newcommand{\angles}[1]{\left\langle #1 \right\rangle}%
\newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}%
\newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}%
\newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}%
\newcommand{\dd}{{\rm d}}%
\newcommand{\ds}[1]{\displaystyle{#1}}%
\newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}%
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}%
\newcommand{\fermi}{\,{\rm f}}%
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}%
\newcommand{\half}{{1 \over 2}}%
\newcommand{\ic}{{\rm i}}%
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}%
\newcommand{\isdiv}{\,\left.\right\vert\,}%
\newcommand{\ket}[1]{\left\vert #1\right\rangle}%
\newcommand{\ol}[1]{\overline{#1}}%
\newcommand{\pars}[1]{\left( #1 \right)}%
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}%
\newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}%
\newcommand{\sech}{\,{\rm sech}}%
\newcommand{\sgn}{\,{\rm sgn}}%
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}%
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
$\ds{I \equiv \int_{0}^{\infty}x^{2}\expo{-x^{2}}{\rm erf}\pars{x}\ln\pars{x}
\,\dd x}$. Let's
$\ds{{\cal I}\pars{\mu}
\equiv \int_{0}^{\infty}x^{\mu}\expo{-x^{2}}{\rm erf}\pars{x}\,\dd x}$ such that $\ds{I = \lim_{\mu \to 2}\totald{{\cal I}\pars{\mu}}{\mu}}$
Since
$\ds{{\rm erf}\pars{x}
\stackrel{{\rm def.}}{=}{2 \over \root{\pi}}\int_{0}^{x}\expo{-y^{2}}\,\dd y}$:
\begin{align}
{\cal I}\pars{\mu}&=\int_{0}^{\infty}x^{\mu}\expo{-x^{2}}
{2 \over \root{\pi}}\int_{0}^{\infty}\Theta\pars{x - y}\expo{-y^{2}}\,\dd y\,\dd x
\\[3mm]&=
{2 \over \root{\pi}}\int_{0}^{\pi/2}\dd\theta\,\cos^{\mu}\pars{\theta}
\Theta\pars{\cos\pars{\theta} - \sin\pars{\theta}}
\overbrace{\int_{0}^{\infty}\dd r\,r^{\mu + 1}\expo{-r^{2}}}
^{\Gamma\pars{1 + \mu/2}/2}
\\[3mm]&={1 \over \root{\pi}}\,\Gamma\pars{1 + {\mu \over 2}}
\int_{0}^{\pi/4}\dd\theta\,\cos^{\mu}\pars{\theta}
\end{align}
\begin{align}
I&=\lim_{\mu \to 2}\totald{{\cal I}\pars{\mu}}{\mu}=
{1 \over 2\root{\pi}}\,\overbrace{\Psi\pars{2}}^{\ds{1 - \gamma}}\
\overbrace{\int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}}^{\ds{\pars{\pi + 2}/8}}
\\[3mm]&+
{1 \over \root{\pi}}\,\overbrace{\Gamma\pars{2}}^{\ds{1}}
\overbrace{%
\int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}\ln\pars{\cos\pars{\theta}}}
^{\ds{\braces{4G + \pi - 2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}}/16}}
\end{align}
$\Gamma\pars{z}$ and $\Psi\pars{z}$ are the $\it Gamma$ and $\it Digamma$ functions, respectively. $\gamma$ and $G$ are the $\it Euler-Mascheroni$ and Catalan constants, respectively.
$$
\begin{array}{|l|}\hline \mbox{}\\
\quad{\displaystyle\int_{0}^{\infty}x^{2}\expo{-x^{2}}
\,\mathrm{erf}\pars{x}\ln\pars{x}\,\dd x}\quad
\\[2mm] =
\quad{{\displaystyle\quad\pars{\pi + 2}\pars{1 - \gamma} + 4G + \pi -
2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}\quad} \over
{\displaystyle 16\root{\pi}}}\quad
\\ \mbox{}\\ \hline
\end{array}
\approx 0.0436462
$$
Best Answer
Consider the integral $$I = \int_0^{\infty} \frac{dx}{(f+c\cdot x)^{a}(1+d\cdot x)^b}$$ where $a,b,c,d,f$ are positive constants as written in the problem.Before I begin, denote $(a)_n = \frac{\Gamma(n+a)}{\Gamma(a)}$ as the Pochhammer symbol of $a$ where $\Gamma$ is the Gamma function, and $$ {}_2F_{1}(a,b,c,z) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n} \frac{z^n}{n!}.$$ Additionally, $\frac{1}{n+z} = \frac{(z)_n}{z(z+1)_n}$ and $$ \left(1+\alpha t\right)^{-a} = \sum_{n=0}^{\infty} \frac{(a)_n}{n!}\left(-\alpha t\right)^n. $$ This follows from the Binomial expansion. Now do the change of variables $1+d\cdot x= \frac{1}{t}$ in I, for $x=0 \mapsto t=1$ and $x\to \infty$ implies $t \to 0$. This results in the Jacobian $-1/(d\cdot t^2)$, and the integral now can be simplified as $$\begin{align*} I &= \frac{1}{d}\int_{0}^{1} \frac{t^{b-2}}{\left(f+\frac{c}{d}\left(\frac{1-t}{t}\right)\right)^a}dt = d^{a-1}\int_{0}^{1} \frac{t^{b+a-2}}{\left(c+(fd-c)t\right)^a}dt,\\ &= \frac{d^{a-1}}{c^{a}}\int_{0}^{1} t^{b+a-2}\left(1+\frac{fd-c}{c}t\right)^{-a}dt \end{align*}$$ Using the Binomial expansion written above, $$\begin{align*} I &= \frac{d^{a-1}}{c^{a}}\int_{0}^{1} \sum_{n=0}^{\infty} \frac{(a)_n}{n!}\left(1-\frac{fd}{c}\right)^{n}t^{n+a+b-2}dt \\ &=\frac{d^{a-1}}{(a+b-1)c^{a}}\sum_{n=0}^{\infty} \frac{(a)_n(a+b-1)_n}{(a+b)_n}\frac{1}{n!}\left(1-\frac{fd}{c}\right)^{n}\end{align*}$$ Thus, $$\boxed{I = \frac{d^{a-1}}{(a+b-1)c^{a}}{}_2F_1\left(a,a+b-1;a+b;1-\frac{fd}{c}\right)}$$
Edit 1: This might seem different from the solution written in the comment for $f=1$. But using Pfaff transformation (see here), the solutions are same since $$ { }_2 F_1(a, b ; c ; z)=(1-z)^{-b}{ }_2 F_1\left(c-a, b; c ; \frac{z}{z-1}\right).$$