On the general form of the family $\sum_{n=1}^\infty \frac{n^{k}}{e^{2n\pi}-1} $ I. $k=4n+3.\;$ From this post, one knows that $$\sum_{n=1}^\infty \frac{n^{3}}{e^{2n\pi}-1} = \frac{\Gamma\big(\tfrac{1}{4}\big)^8}{2^{10}\cdot5\,\pi^6}-\frac{1}{240}$$ and a Mathematica session reveals$$\sum_{n=1}^\infty  \frac{n^{7}}{e^{2n\pi}-1} =\frac{3\,\Gamma\big(\tfrac{1}{4}\big)^{16}}{2^{17}\cdot5\,\pi^{12}}-\frac{1}{480}$$
$$\sum_{n=1}^\infty  \frac{n^{11}}{e^{2n\pi}-1} =\frac{189\,\Gamma\big(\tfrac{1}{4}\big)^{24}}{2^{22}\cdot5\cdot13\,\pi^{18}}-\frac{691}{65520}$$
and so on. The $691$ is a clue that Bernoulli numbers are involved. 
II. $k=4n+1.\;$ It evaluates to a rational number,
$$\sum_{n=1}^\infty  \frac{n^{5}}{e^{2n\pi}-1} =\frac{1}{504}$$
$$\sum_{n=1}^\infty  \frac{n^{9}}{e^{2n\pi}-1} =\frac{1}{264}$$
$$\sum_{n=1}^\infty  \frac{n^{13}}{e^{2n\pi}-1} =\frac{1}{24}$$ 
etc, with the last mentioned in this post.
Q: What are the general forms of I and II in terms of the Bernoulli numbers? (And a reference to Ramanujan's Notebooks, if possible.)


Edit added June 2019.

Since the gamma function/pi ratio involved has a simple form in terms of the elliptic integral singular value $K(k_n)$,
$$\beta_1=\frac{\Gamma\big(\tfrac14\big)^8}{2^8\pi^6}=\left(\frac{K(k_1)}{\pi}\right)^4$$
hence was the case $\tau=\sqrt{-1}$. So by analogy, the case  $\tau=\sqrt{-3}$,
$$\beta_3=\frac{3\Gamma\big(\tfrac13\big)^{12}}{2^9\cdot2^{1/3}\,\pi^8}=\left(\frac{K(k_3)}{\pi}\right)^4$$
And a quick test showed,
$$\sum_{k=1}^\infty \frac{k^{3}}{e^{2\pi\sqrt3\,k}-1} = \frac{1}{16}\color{red}{\beta_3}-\frac{1}{240}$$ 
$$\sum_{k=1}^\infty \frac{k^{7}}{e^{2\pi\sqrt3\,k}-1} = \frac{17}{32}\color{red}{{\beta_3}^2}-\frac{1}{480}$$
and so on. 
 A: You are asking about the special value $E_{2k}(i)$ of Eisenstein series $E_{2k}(z)$ and more generally $f(i)$ for any modular form $\in  M_k(SL_2(\Bbb{Z}))$.

Any modular form $\in M_k(SL_2(\Bbb{Z}))$ is in $\Bbb{C}[E_4,E_6]$.

Both $\Delta(z) = q \prod_{n \ge 1}(1-q^n)^{24}$ and $E_4(z)^3-E_6(z)^2$ are non-zero cusp forms $\in S_{12}(SL_2(\Bbb{Z})$, from $\dim(S_{12}(SL_2(\Bbb{Z}))=1$ we obtain that $E_4(z)^3-E_6(z)^2 = c\Delta(z)$. 
$\Delta(z)$ has only one simple zero at $i\infty$ so $f \mapsto f \Delta$ is an isomorphism $M_k(SL_2(\Bbb{Z}))\to S_{k+12}(SL_2(\Bbb{Z})$ so that $M_{k+12}(SL_2(\Bbb{Z})= \Bbb{C} E_4^a E_6^b+ \Delta M_k(SL_2(\Bbb{Z})$ and it suffices to check the claim for $k < 12$ where it reduces to that $\dim(M_k(SL_2(\Bbb{Z})) \le 1$.

For $k$ odd then $E_{2k}(i) = i^{-2k} E_{2k}(-1/i) = -E_{2k}(-1/i)=0$. 

$ $

The elliptic integrals corresponding to $\Bbb{Z}+i\Bbb{Z}$ reduce to evaluating the beta function from which we find $\eta(i) = \frac{\Gamma(1/4)}{2 \pi^{3/4}}$ (can someone help to fill that claim ?) and $E_4(i)= 240 \frac{\Gamma\left(\frac{1}{4}\right)^8}{2^{10}\cdot 5\,\pi^6}$. 

Thus it reduces to find the polynomial $P(x,y) \in \Bbb{Q}[x,y]$ such that $E_{2k}(z) = P(E_4(z),E_6(z))$. 
$\Bbb{C}[E_4,E_6]$ is graded by weight of modularity so that $$E_{2k}(z) = \sum_{4a+6b=2k} c_a E_4(z)^aE_6(z)^b$$  the $c_a \in \Bbb{Q}$ are found from solving a linear equation in the first few $q$-expansion coefficients of $E_{2k}(z),E_4(z)^aE_6(z)^b$
$$E_{2k}(i) = P(E_4(i),E_6(i))= \cases{0 \text{ if } 2  \nmid k\\ c_{k/2} \left(240\frac{\Gamma\left(\frac{1}{4}\right)^8}{2^{10}\cdot 5\,\pi^6}\right)^{k/2} \text{ otherwise}}$$

There is a possible algorithm : say $2k = 12M+4a+6b, 4a+6b < 12$ and $$g_M(z) = E_{2k}(z), \qquad g_{m-1}(z) = \frac{1}{\Delta(z)} (g_m(z)-g_m(i\infty) E_4(z)^a E_6(z)^{b+2m}) $$ Then $g_m \in M_{12m+l}(SL_2(\Bbb{Z}))$,
 $g_0(z) = g_0(i\infty) E_4(z)^aE_6(z)^b$, $E_{2k}(z) = \sum_{m=0}^M \Delta(z)^{M-m} g_m(i\infty) E_4(z)^a E_6(z)^{b+2m}$
$$E_{2k}(i) =\Delta(i)^M g_0(i\infty)= E_4(i)^{3M} g_0(i\infty), \qquad g_0(i\infty) \in \frac{1}{d_k}\Bbb{Z}$$
where $d_k$ is the numerator of $B_{2k}$.
A: Ramanujan deals with the sums of type $$S_{r}(q) = -\frac{B_{r + 1}}{2(r + 1)} + \sum_{n = 1}^{\infty}\frac{n^{r}q^{n}}{1 - q^{n}}\tag{1}$$ where $B_{r}$ are Bernoulli's Numbers defined by $$\frac{x}{e^{x} - 1} = \sum_{r = 0}^{\infty}B_{r}\frac{x^{r}}{r!}\tag{2}$$ in his paper On certain arithmetical functions which appeared in Transactions of the Cambridge Philosophical Society in 1916. And he gives the general recursion formula for $S_{r}$ through which we can calculate the value of $S_{r}(q)$ in terms of a polynomial in functions $P(q), Q(q), R(q)$ which are given by $$P(q) = -24S_{1}(q), Q(q) = 240S_{3}(q), R(q) = -504S_{5}(q)\tag{3}$$ Moreover using link between theta functions and elliptic integrals it is possible to express $P, Q, R$ in terms of elliptic integral $K$ and modulus $k$ where $k$ corresponds to nome $q$. It is thus possible to express $S_{r}(q)$ as a polynomial in $K, k$. The sum in your question deals with the values of $$S_{r}(q^{2}) + \frac{B_{r + 1}}{2(r + 1)}$$ for $q = e^{-\pi}$ which translates to $k = 1/\sqrt{2}$ and $K = \Gamma^{2}(1/4)/4\sqrt{\pi}$. And hence one should expect the occurrence of $\Gamma(1/4)$ in the evaluations. When $r$ is of type $4m + 1$ then the value of $S_{r}$ is always a polynomial with factor $R$ (this is proved by Ramanujan) and for $q = e^{-\pi}$ the value of $R(q^{2})$ is $0$ because $$R(q^{2}) = \left(\frac{2K}{\pi}\right)^{6}(1 + k^{2})(1 - 2k^{2})\left(1 - \frac{k^{2}}{2}\right)$$ which vanishes when $k = 1/\sqrt{2}$. It follows that for $r = 4m + 1$ we have the desired sum as $\dfrac{B_{r + 1}}{2(r + 1)}$ (this gives the general formula for sum of type II in your question).
When $r$ is of type $r = 4m + 3$ the desired sum is expressed as a rational number plus some expression consisting of $\Gamma(1/4)$ and $\pi$. Again the rational number in this expression is $B_{r + 1}/2(r + 1)$. A general formula is not known, but using Ramanujan's table of values of $S_{r}(q)$ in his paper we can get this for all odd $r$ upto $r = 31$.
Thus for example Ramanujan gives the formula for $r = 31$ as $$7709321041217 + 32640\sum_{n = 1}^{\infty}\frac{n^{31}q^{n}}{1 - q^{n}} = 764412173217Q^{8}(q) + \text{ terms containing }R(q)\tag{4}$$ and therefore $$\sum_{n = 1}^{\infty}\frac{n^{31}}{e^{2\pi n} - 1} = \frac{764412173217Q^{8}(e^{-2\pi})}{32640} - \frac{7709321041217}{32640}$$ where $$Q(e^{-2\pi}) = Q(q^{2}) = \left(\frac{2K}{\pi}\right)^{4}(1 - k^{2} + k^{4})$$ with $k = 1/\sqrt{2}$ and $K = \Gamma^{2}(1/4)/4\sqrt{\pi}$.
Also note that the expression for $S_{r}$ for $r = 4m + 3$ has only one term without $R$ and that is a rational multiple of $Q^{(r + 1)/4}$ and therefore it follows that the sum in question for $r = 4m + 3$ is of the form $$A\cdot\frac{\Gamma^{2(r + 1)}(1/4)}{\pi^{3(r + 1)/2}} + \frac{B_{r + 1}}{2(r + 1)}$$ where $A$ is a rational number.

An exposition of Ramanujan's paper mentioned above is given in my blog posts here and here.
