Suppose we seek to evaluate
$$F(p) = \sum_{n\ge 1} \frac{n^{4p-1}}{e^{\pi n}-1}
- 16^p \sum_{n\ge 1} \frac{n^{4p-1}}{e^{4\pi n}-1}.$$
These sums may be evaluated
using harmonic summation techniques.
Introduce the sum
$$S(x; p) = \sum_{n\ge 1} \frac{n^{4p-1}}{e^{nx}-1}$$
with $p$ a positive integer and $x\ge 0.$
The sum term is harmonic and may be evaluated by inverting its Mellin
transform.
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = k^{4p-1}, \quad \mu_k = k
\quad \text{and} \quad
g(x) = \frac{1}{e^x-1}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$\int_0^\infty \frac{1}{e^{x}-1} x^{s-1} dx
= \int_0^\infty \frac{e^{-x}}{1-e^{-x}} x^{s-1} dx
\\ = \int_0^\infty \sum_{q\ge 1} e^{-q x} x^{s-1} dx
= \sum_{q\ge 1} \int_0^\infty e^{-q x} x^{s-1} dx
\\= \Gamma(s) \sum_{q\ge 1} \frac{1}{q^s}
= \Gamma(s) \zeta(s).$$
It follows that the Mellin transform $Q(s)$ of the harmonic sum
$S(x,p)$ is given by
$$Q(s) = \Gamma(s) \zeta(s) \zeta(s-(4p-1))
\\ \text{because}\quad
\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} =
\sum_{k\ge 1} k^{4p-1} \frac{1}{k^s}
= \zeta(s-(4p-1))$$
for $\Re(s) > 4p.$
The Mellin inversion integral here is
$$\frac{1}{2\pi i} \int_{4p+1/2-i\infty}^{4p+1/2+i\infty} Q(s)/x^s ds$$
which we evaluate by shifting it to the left for an expansion about
zero.
The two zeta function terms cancel the poles of the gamma function
term and we are left with just
$$\begin{align}
\mathrm{Res}(Q(s)/x^s; s=4p) & = \Gamma(4p) \zeta(4p)
\quad\text{and}\\
\mathrm{Res}(Q(s)/x^s; s=0) & = \zeta(0) \zeta(-(4p-1)).
\end{align}$$
Computing these residues we get
$$- (4p-1)! \frac{B_{4p} (2\pi)^{4p}}{2(4p)!}
= - \frac{B_{4p} (2\pi)^{4p}}{2\times 4p}
\quad\text{and}\quad
- \frac{1}{2} \times -\frac{B_{4p}}{4p}.$$
This shows that
$$S(x;p) = - \frac{B_{4p} (2\pi)^{4p}}{8p\times x^{4p}}
+ \frac{B_{4p}}{8p}
+ \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} Q(s)/x^s ds.$$
To treat the integral recall the duplication formula of the gamma
function:
$$\Gamma(s) =
\frac{1}{\sqrt\pi} 2^{s-1}
\Gamma\left(\frac{s}{2}\right)
\Gamma\left(\frac{s+1}{2}\right).$$
which yields for $Q(s)$
$$\frac{1}{\sqrt\pi} 2^{s-1}
\Gamma\left(\frac{s}{2}\right)
\Gamma\left(\frac{s+1}{2}\right)
\zeta(s) \zeta(s-(4p-1))$$
Furthermore observe the following variant of the functional equation
of the Riemann zeta function:
$$\Gamma\left(\frac{s}{2}\right)\zeta(s)
= \pi^{s-1/2} \Gamma\left(\frac{1-s}{2}\right)
\zeta(1-s)$$
which gives for $Q(s)$
$$\frac{1}{\sqrt\pi} 2^{s-1}
\pi^{s-1/2}
\Gamma\left(\frac{s+1}{2}\right)
\Gamma\left(\frac{1-s}{2}\right)
\zeta(1-s)\zeta(s-(4p-1))
\\ =
\frac{1}{\sqrt\pi} 2^{s-1}
\pi^{s-1/2}
\frac{\pi}{\sin(\pi(s+1)/2)}
\zeta(1-s)\zeta(s-(4p-1))
\\ =
2^{s-1}
\frac{\pi^s}{\sin(\pi(s+1)/2)}
\zeta(1-s)\zeta(s-(4p-1)).$$
Now put $s=4p-u$ in the remainder integral to get
$$\frac{1}{x^{4p}}
\frac{1}{2\pi i} \int_{4p+1/2-i\infty}^{4p+1/2+i\infty}
2^{4p-1-u}
\\ \times \frac{\pi^{4p-u}}{\sin(\pi(4p+1-u)/2)}
\zeta(u-(4p-1))\zeta(1-u) x^u du
\\ = \frac{2^{4p} \pi^{4p}}{x^{4p}}
\frac{1}{2\pi i} \int_{4p+1/2-i\infty}^{4p+1/2+i\infty}
2^{u-1}
\\ \times \frac{\pi^{u}}{\sin(\pi(4p+1-u)/2)}
\zeta(u-(4p-1))\zeta(1-u) (x/\pi^2/2^2)^u du.$$
Now $$\sin(\pi(4p+1-u)/2) = \sin(\pi(1-u)/2+2\pi p)
\\ = \sin(\pi(1-u)/2) = - \sin(\pi(-1-u)/2)
= \sin(\pi(u+1)/2).$$
We have shown that
$$\bbox[5px,border:2px solid #00A000]
{S(x;p) = - \frac{B_{4p} (2\pi)^{4p}}{8p\times x^{4p}}
+ \frac{B_{4p}}{8p}
+ \frac{2^{4p} \pi^{4p}}{x^{4p}} S(4\pi^2/x;p)}.$$
In particular we get
$$S(\pi;p) = - \frac{B_{4p} 2^{4p}}{8p}
+ \frac{B_{4p}}{8p}
+ 2^{4p} S(4\pi;p)$$
and
$$S(4\pi;p) = - \frac{B_{4p} 2^{-4p}}{8p}
+ \frac{B_{4p}}{8p}
+ 2^{-4p} S(\pi;p).$$
Therefore
$$S(\pi;p)-2^{4p}S(4\pi;p) \\ =
- \frac{B_{4p} (2^{4p}-1)}{8p}
+ \frac{B_{4p}}{8p} (1-2^{4p})
+ (2^{4p} S(4\pi;p)-S(\pi;p)).$$
We finally conclude that
$$F(p) = \frac{B_{4p}}{4p} (1-2^{4p}) - F(p)$$
or
$${\large\color{green}{F(p) =
\frac{B_{4p}}{8p} (1-2^{4p})}}.$$
The first few values are
$$\frac{1}{16},{\frac {17}{32}},{\frac {691}{16}},
{\frac {929569}{64}},{\frac {221930581}{16}},\ldots$$
We can use the functional equation to extract additional sum formulae.
For example when we choose the pair
$$\sqrt{2}\pi \quad\text{and}\quad 2\sqrt{2}\pi$$
the scalar is $2^{2p}.$ The calculation starts from
$$S(\sqrt{2}\pi;p) = - \frac{B_{4p} 2^{2p}}{8p}
+ \frac{B_{4p}}{8p}
+ 2^{2p} S(2\sqrt{2}\pi;p)$$
and
$$S(2\sqrt{2}\pi;p) = - \frac{B_{4p} 2^{-2p}}{8p}
+ \frac{B_{4p}}{8p}
+ 2^{-2p} S(\sqrt{2}\pi;p).$$
Therefore
$$S(\sqrt{2}\pi;p)-2^{2p}S(2\sqrt{2}\pi;p) \\ =
- \frac{B_{4p} (2^{2p}-1)}{8p}
+ \frac{B_{4p}}{8p} (1-2^{2p})
+ (2^{2p} S(2\sqrt{2}\pi;p)-S(\sqrt{2}\pi;p)).$$
We obtain the formula
$$\sum_{n\ge 1} \frac{n^{4p-1}}{e^{\sqrt{2}\pi n}-1}
- 4^p \sum_{n\ge 1} \frac{n^{4p-1}}{e^{2\sqrt{2}\pi n}-1}
= \frac{B_{4p}}{8p} (1-2^{2p}).$$
We get for the initial values (factor is $1+2^{2p}$)
$${\frac {1}{80}},\frac{1}{32},{\frac {691}{1040}},
{\frac {3617}{64}},{\frac {5412941}{400}},\ldots$$
Addendum. To illustrate the creation of identities from the
functional equation we present a third example, which is
$$\sqrt{3}\pi\quad\text{and}\quad 4\pi/\sqrt{3}.$$
In this example the scalar is $2^{4p} 3^{-2p}.$
The calculation starts from
$$S(\sqrt{3}\pi;p) = - \frac{B_{4p} 2^{4p}}{8p\times 3^{2p}}
+ \frac{B_{4p}}{8p}
+ \frac{2^{4p}}{3^{2p}} S(4\pi/\sqrt{3};p)$$
and
$$S(4\pi/\sqrt{3};p) = - \frac{B_{4p} 3^{2p}}{8p\times 2^{4p}}
+ \frac{B_{4p}}{8p}
+ \frac{3^{2p}}{2^{4p}} S(\sqrt{3}\pi;p).$$
Therefore
$$S(\sqrt{3}\pi;p)-2^{4p} 3^{-2p} S(4\pi/\sqrt{3};p) \\ =
- \frac{B_{4p} (2^{4p} 3^{-2p}-1)}{8p}
+ \frac{B_{4p}}{8p} (1-2^{4p} 3^{-2p})
+ (2^{4p} 3^{-2p} S(\sqrt{3}\pi;p)-S(4\pi/\sqrt{3};p)).$$
We obtain the formula
$$\sum_{n\ge 1} \frac{n^{4p-1}}{e^{n \sqrt{3}\pi}-1}
- 2^{4p} 3^{-2p}
\sum_{n\ge 1} \frac{n^{4p-1}}{e^{n 4\pi/\sqrt{3}}-1}
= \frac{B_{4p}}{8p} (1-2^{4p} 3^{-2p}).$$
The reader is cordially invited to prove this last result by a
different method.
Remark. The general pattern for
$$\beta\quad\text{and}\quad 4\pi^2/\beta$$
is
$$\bbox[5px,border:2px solid #00A000]
{\sum_{n\ge 1} \frac{n^{4p-1}}{e^{n \beta}-1}
- \frac{2^{4p}\pi^{4p}}{\beta^{4p}}
\sum_{n\ge 1} \frac{n^{4p-1}}{e^{n 4\pi^2/\beta}-1}
= \frac{B_{4p}}{8p}
\left(1-\frac{2^{4p}\pi^{4p}}{\beta^{4p}}\right).}$$
