The closed-form solution of the family $\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(pn+m)}$?
$$\boxed{\int_0^1 \frac{\ln(1-x) \ln(1-x^5)}{x} dx=\\ 4\zeta(3)-\frac{\pi}{5}\operatorname{Cl}_2\left(\frac{4\pi}{5}\right)-\frac{3\pi}{5}\operatorname{Cl}_2\left(\frac{2\pi}{5}\right)+3\operatorname{Cl}_3\left(\frac{4\pi}{5}\right)+3\operatorname{Cl}_3\left(\frac{2\pi}{5}\right)}$$ $$\operatorname{Cl}_2\left(x\right)=\sum_{n=1}^\infty \frac{\sin(nx)}{n^2},\quad \operatorname{Cl}_3\left(x\right)=\sum_{n=1}^\infty \frac{\cos(nx)}{n^3}$$
(Added by OP.) But since,
$$\operatorname{Cl}_3\left(\frac{4\pi}{5}\right)+\operatorname{Cl}_3\left(\frac{2\pi}{5}\right) =-\frac{12}{25}\zeta(3)$$
then the above can be simplified as,
$$\boxed{\int_0^1 \frac{\ln(1-x) \ln(1-x^5)}{x} dx=\frac{64}{25}\zeta(3)-\frac{\pi}{5}\operatorname{Cl}_2\left(\frac{4\pi}{5}\right)-\frac{3\pi}{5}\operatorname{Cl}_2\left(\frac{2\pi}{5}\right)}$$
Tools used: $$(1-x^5)=(1-x)(1+\varphi x+x^2)(1-\frac{1}{\varphi}x+x^2), \quad \varphi =\frac{\sqrt 5+1}{2} $$ $$\ln(1+\varphi x+x^2)=-2\sum_{n=1}^\infty \frac{\cos\left(\frac{4n\pi}{5}\right)}{n}x^n$$ $$\ln(1-\frac{1}{\varphi} x+x^2)=-2\sum_{n=1}^\infty \frac{\cos\left(\frac{2n\pi}{5}\right)}{n}x^n$$ $$\int_0^1 x^{n-1}\ln(1-x)dx=-\frac1n\sum_{k=1}^n \frac{1}{k}=-\frac{H_n}{n}$$ $$S(x)=\sum_{n=1}^\infty \frac{x^n}{n^2}H_n=\operatorname{Li}_3(x)-\operatorname{Li}_3(1-x)+\operatorname{Li}_2(1-x)\ln(1-x)+\frac{1}{2}\ln x \ln^2(1-x)+\zeta(3) $$
$$\small I(5)=\int_0^1 \frac{\ln^2(1-x)}{x}dx+\int_0^1 \frac{\ln(1-x)\ln(1+\varphi x+x^2)}{x}dx+\int_0^1\frac{\ln(1-x)\ln(1-\frac{1}{\varphi} x+x^2)}{x}dx$$ $$=\sum_{n=1}^\infty \int_0^1 x^{n-1} \ln^2 xdx-2\sum_{n=1}^\infty \frac{\cos\left(\frac{4n\pi}{5}\right)+\cos\left(\frac{2n\pi}{5}\right)}{n}\int_0^1 x^{n-1} \ln(1-x)dx$$ $$=2\sum_{n=1}^\infty \frac{1}{n^3}+2\sum_{n=1}^\infty \frac{\cos\left(\frac{4n\pi}{5}\right)+\cos\left(\frac{2n\pi}{5}\right)}{n^2}H_n=2\zeta(3)+2\Re \left(S\left(e^{4 i \pi/5}\right)+S\left(e^{2 i \pi/5}\right)\right)\tag 1$$
In order to calculate the real parts of the polylogs I used this approach to find: $$\Re \operatorname{Li}_3(e^{4i\pi/5})=\operatorname{Cl}_3\left(\frac{4\pi}{5}\right)$$ $$\Re \operatorname{Li}_3(1-e^{4i\pi/5})=\frac{\zeta(3)}{2}-\frac12 \operatorname{Cl}_3\left(\frac{4\pi}{5}\right)+\frac{2\pi^2}{25}\ln\left(\frac{5+\sqrt 5}{2}\right)$$ $$\Re \operatorname{Li}_2(1-e^{4i\pi/5})\ln(1-e^{i4\pi/5})=\frac{\pi^2}{25}\ln\left(\frac{5+\sqrt 5}{2}\right)-\frac{\pi}{10}\operatorname{Cl}_2\left(\frac{4\pi}{5}\right)$$ $$\Re \ln(e^{i4\pi/5})\ln^2(1-e^{i4\pi/5})=\frac{2\pi^2}{25}\ln\left(\frac{5+\sqrt 5}{2}\right)$$
$$\Re \operatorname{Li}_3(e^{2i\pi/5})=\operatorname{Cl}_3\left(\frac{2\pi}{5}\right)$$ $$\Re \operatorname{Li}_3(1-e^{2i\pi/5})=\frac{\zeta(3)}{2}-\frac12 \operatorname{Cl}_3\left(\frac{2\pi}{5}\right)+\frac{\pi^2}{50}\ln\left(\frac{5-\sqrt 5}{2}\right)$$ $$\Re \operatorname{Li}_2(1-e^{4i\pi/5})\ln(1-e^{i4\pi/5})=-\frac{\pi^2}{25}\ln\left(\frac{5-\sqrt 5}{2}\right)-\frac{3\pi}{10}\operatorname{Cl}_2\left(\frac{2\pi}{5}\right)$$ $$\Re \ln(e^{i4\pi/5})\ln^2(1-e^{i4\pi/5})=\frac{3\pi^2}{25}\ln\left(\frac{5-\sqrt 5}{2}\right)$$
And plugging those values in $(1)$ yields the announced result.
We may apply a discrete Fourier transform to the following generating function $$\sum_{n=1}^\infty \frac{x^n}{n^2}H_n=\operatorname{Li}_3(x)-\operatorname{Li}_3(1-x)+\operatorname{Li}_2(1-x)\ln(1-x)+\frac{1}{2}\ln x \ln^2(1-x)+\zeta(3)$$ since $$ I(p) = \sum_{n\geq 1}\frac{H_{p n}}{pn^2}. $$ The only term leading to a non-elementary contribution is the sum of $\operatorname{Li}_3(1-x)$ over the $p$-th roots of unity.
Some work in progress on the general series. No closed form, sorry, but I think this may be interesting anyway.
Let's study the function $I(p)$. Obviously:
$$I \left( \frac{1}{p} \right)= p I(p)$$
Thus we are only interested in the case $p \geq 1$.
Let's sum over $m$. This gives us:
$$I(p)=\frac{\pi^2}{6}\frac{\gamma}{p}+\frac{1}{p} \sum_{n=1}^\infty \frac{\psi(pn+1)}{n^2} \tag{1}$$
There's a lot of various identities for polygamma which could be useful here.
1) Consider the following identity:
$$\psi(pn+1)=\log (pn+1)-\sum_{k=1}^\infty \frac{|G_k| (k-1)!}{(pn+1)_k}$$
Where $G_k$ are so called Gregory coefficients. $G_k= \int_0^1 \binom{x}{k} dx$ and $|G_k| \asymp \frac{1}{k \log^2 k}$ if $k \to \infty$.
$$I(p)=\frac{\pi^2}{6}\frac{\gamma+\log p}{p}+\frac{1}{p} \sum_{n=1}^\infty \frac{\log(n+1/p)}{n^2}-\frac{1}{p} \sum_{k=1}^\infty \frac{|G_k| k!}{k} \sum_{n=0}^\infty \frac{1}{(n+1)^2 (pn+p+1)_k} $$
The second series doesn't have a closed form as far as I know but it's elementary at least.
Third double series should be small in value and you might notice I changed the order of summation.
$$\sum_{n=0}^\infty \frac{1}{(n+1)^2 (pn+p+1)_k}= \frac{p!}{(p+k)!} {_{k+3} F_{k+2}} \left( \begin{array}(1,1,1, \frac{1}{p}+1, \ldots, \frac{1}{p}+k \\ 2,2,\frac{1}{p}+2, \ldots, \frac{1}{p}+k+1 \end{array};1 \right)$$
So we have:
$$pI(p)=\frac{\pi^2}{6}(\gamma+\log p)+\sum_{n=1}^\infty \frac{\log(n+1/p)}{n^2}- \\ -\sum_{k=1}^\infty \frac{|G_k|}{k \binom{p+k}{k}} {_{k+3} F_{k+2}} \left( \begin{array}(1,1,1, \frac{1}{p}+1, \ldots, \frac{1}{p}+k \\ 2,2,\frac{1}{p}+2, \ldots, \frac{1}{p}+k+1 \end{array};1 \right) \tag{2}$$
For $p>1$ the first terms and the log series give the most important contribution. The last series is complicated, but we can easily compute any finite number of terms to get more digits.
Further expanding the logarithm and using:
$$\sum_{n=1}^\infty \frac{\log(n)}{n^2}=- \frac{\pi^2}{6} (\gamma+ \log(2 \pi))+2 \pi^2 \log A $$
Whee A is http://mathworld.wolfram.com/Glaisher-KinkelinConstant.html.
We have:
$$pI(p)=\frac{\pi^2}{6}(\log p+12 \log A-\log 2 \pi)+\sum_{n=1}^\infty \frac{1}{n^2} \log \left(1+\frac{1}{pn} \right)- \\ -\sum_{k=1}^\infty \frac{|G_k|}{k \binom{p+k}{k}} {_{k+3} F_{k+2}} \left( \begin{array}(1,1,1, \frac{1}{p}+1, \ldots, \frac{1}{p}+k \\ 2,2,\frac{1}{p}+2, \ldots, \frac{1}{p}+k+1 \end{array};1 \right) \tag{3}$$
For $p \to \infty$ the asymptotic expansion will then be:
$$p I(p) \asymp \frac{\pi^2}{6}(\log p+12 \log A-\log 2 \pi)+ \frac{\zeta(3)}{2p} \tag{4}$$
Where an additional $-\zeta(3)/(2p)$ comes from the third series as the first term in asymptotic expansion for large $p$.
An exapmle:
$$100 I(100)=9.4682325532367113866$$
$$\frac{\pi^2}{6}(\log 100+12 \log A-\log 2 \pi)+ \frac{\zeta(3)}{2 \cdot 100}=9.4682415725122177074876$$
As you can see the asymptotic expansion works well enough, though some further correction terms are needed.
From (1), expanding the logarithm as we did, and using the well known asymptotic expansion of harmonic numbers we can make a full asymptotic series:
$$p I(p) \asymp \frac{\pi^2}{6}(\log p+12 \log A-\log 2 \pi)+\frac{\zeta(3)}{2p} -\sum_{k=1}^\infty \frac{B_{2k}}{2k p^{2k}} \zeta(2k+2) \tag{5}$$
I'll check it numerically later, but I'm pretty sure it doesn't converge. Still, for large $p$ a few first terms should give a lot of correct digits.
Using the explicit form for even zetas, we have:
$$p I(p) \asymp \frac{\pi^2}{6}\log \frac{p}{2\pi}+2\pi^2 \log A+\frac{\zeta(3)}{2p} -\frac{\pi^2}{2} \sum_{k=1}^\infty \frac{(-1)^k B_{2k}B_{2k+2}}{k(k+1) (2k+1)!} \frac{(2\pi)^{2k}}{p^{2k}} \tag{6}$$
The logarithmic terms and the series make me think that $p=2\pi$ is some special value.