Solution (source code)

= Solution

Write $\boldsymbol\sigma_i=\mathbf m+\delta\boldsymbol\sigma_i$ and neglect products of two fluctuations. Then
$$
\boldsymbol\sigma_i\cdot\boldsymbol\sigma_j
\simeq\mathbf m\cdot(\boldsymbol\sigma_i+\boldsymbol\sigma_j)-m^2.
$$
Each site has $q$ neighbours and each bond is counted once, so the <mean-field approximation> gives
$$
E_{\rm MF}=-Jq\,\mathbf m\cdot\sum_i\boldsymbol\sigma_i
+\frac{NJq}{2}m^2.
$$
Choose the $x$ axis along $\mathbf m$. For the <four-state clock model>, the single-site <partition function> is
$$
z_1=\sum_{\theta=0,\pi/2,\pi,3\pi/2}e^{\beta Jqm\cos\theta}
=2+2\cosh(\beta Jqm).
$$
Consequently
$$
f=\frac{Jq}{2}m^2-T\log[2+2\cosh(\beta Jqm)],
$$
and therefore
$$
\boxed{A=\frac{Jq}{2},\qquad C=2,\qquad D=2}.
$$