(1)-(3), proposed for the description of FeBr$_2$), $$\begin{aligned}
{\cal H}^1_{nd}(J^1_{xz}) &=& - J^1_{xz} \sum\limits_{ijk} \left[ S^z_{ijk} S^x_{i+1 \, jk}
- S^x_{ijk} S^z_{i+1 \, jk} \right. 
\nonumber \\
&& \left. + S^z_{ijk} S^y_{i \, j+1 \, k} - S^y_{ijk} S^z_{i \, j+1 \, k}
\right],\end{aligned}$$ or (2) a term involving products of squares of spins, being invariant under the symmetry operations of the space group, $${\cal H}^2_{nd}(J^2_{xz}) = - J^2_{xz} \sum\limits_{\left< \alpha \beta \right>} 
\left( S^z_{\alpha} \right)^2 \left[ \left( S^x_{\beta} \right)^2 + \left( S^y_{\beta} \right)^2
\right],$$ where the sum is over nearest–neighbor in–layer sites.

I consider in the following S = 1 spins where the quantum nature is treated in a semi–classical approximation: the $z$–component is discretized and can only take the values 1, 0, or $-1$, whereas the spin length is fixed to be $\left| \vec{S} \right| = \sqrt{\mbox{S} \left( \mbox{S} + 1 \right)} = \sqrt{2}$. Hence, the spins rotate in the $xy$ plane like a classical vector. The $xy$–components provide additional energy contributions as compared to the S = 1 Ising model, even at low temperatures. Note that in Ref. 9 the spin length was fixed at 1, the spins thus having a planar component only if $S^z = 