Calculation of gain at certain frequency in case of nth order IIR filter
You just need to evaluate the transfer function on the unit circle at the frequency of interest:
$$H(e^{j\omega_0})=\frac{\displaystyle\sum_{k=0}^Nb_ke^{-jk\omega_0}}{\displaystyle\sum_{k=0}^Na_ke^{-jk\omega_0}}\tag{1}$$
and take the magnitude.
For the special values $\omega_0=0$ and $\omega_0=\pi$, Eq. $(1)$ simplifies to
$$H(e^{j0})=\frac{\displaystyle\sum_{k=0}^Nb_k}{\displaystyle\sum_{k=0}^Na_k}\tag{2}$$
and
$$H(e^{j\pi})=\frac{\displaystyle\sum_{k=0}^N(-1)^kb_k}{\displaystyle\sum_{k=0}^N(-1)^ka_k}\tag{3}$$
respectively.
EDIT: If you want a formula that directly expresses the squared magnitude of $H(e^{j\omega})$ then use this:
$$\big|H(e^{j\omega})\big|^2=\frac{\displaystyle r_b[0]+2\sum_{k=1}^Nr_b[k]\cos(k\omega)}{\displaystyle r_a[0]+2\sum_{k=1}^Nr_a[k]\cos(k\omega)}\tag{4}$$
where $r_a[k]$ and $r_b[k]$ are the autocorrelations of the denominator and numerator coefficients, respectively:
$$r_a[k]=a[k]\star a[-k]\\r_b[k]=b[k]\star b[-k]$$
where $\star$ denotes convolution.