Solution (source code)

= Solution

Here $p(n)$ is the <probability> of microstate $n$, including all occupation data. In the <grand canonical ensemble>, with energy $E_n$, particle number $N_n$, <chemical potential> $\mu$ and $T>0$,
$$
p(n)=\mathcal Z^{-1}e^{-(E_n-\mu N_n)/(kT)},\qquad
\mathcal Z=\sum_ne^{-(E_n-\mu N_n)/(kT)}.
$$
Substitution in the <Gibbs entropy> gives $S=k\log\mathcal Z+(U-\mu N)/T$, where $U=\langle E\rangle$ and $N=\langle N_n\rangle$. Differentiating at fixed $\mu,V$ gives $\partial_T\log\mathcal Z=(U-\mu N)/(kT^2)$. Therefore
$$
\boxed{S=k\left[\frac{\partial}{\partial T}(T\log\mathcal Z)\right]_{\mu,V}.}
$$