  Free energy, internal energy, and specific heat for each of the three phonon
spectrum branches of a magnetized Coulomb crystal with body-centered cubic
lattice are calculated by numerical integration over the Brillouin zone in the
range of magnetic fields $B$ and temperatures $T$, such that $0 \le \omega_{\rm
B}/\omega_{\rm p}\le 10^3$ and $10^{-4} \le T/T_{\rm p} \le 10^4$. In this
case, $\omega_{\rm B}$ is the ion cyclotron frequency, $\omega_{\rm p}$ and
$T_{\rm p}$ are the ion plasma frequency and plasma temperature, respectively.
The results of numerical calculations are approximated by simple analytical
formulas. For illustration, these formulas are used to analyze the behavior of
the heat capacity in the crust of a neutron star with strong magnetic field.
Thermodynamic functions of magnetized neutron star crust are needed for
modeling various observational phenomena in magnetars and high magnetic field
pulsars.
