We consider a body centered cubic crystal with one atom in the basis. The primitive lattice vectors are,

\begin{equation} \vec{a}_1=\frac{a}{2}(\hat{x}+\hat{y}-\hat{z}),\quad \vec{a}_2=\frac{a}{2}(-\hat{x}+\hat{y}+\hat{z}),\quad\vec{a}_3=\frac{a}{2}(\hat{x}-\hat{y}+\hat{z}). \end{equation}

Each atom has 8 nearest neighbors and 6 next-nearest neighbors. Let $u^x_{l,m,n}$ the displacement of an atom in the unit cell at position $\vec{r}=l\vec{a}_1+m\vec{a}_2+n\vec{a}_3$ in the $x$-direction from its equilibrium position. If the restoring force that pushes the atoms back to their equilibrium positions is modelled with linear springs with spring constants $C_{1}$ for nearest neighbors and $C_{2}$ for next-nearest neighbors, then the equations of motion for the atoms are,

\begin{equation} M\frac{d^2u_{l,m,n}^{x}}{dt^2}=\frac{C_1}{3}\Big[(u_{l+1,m,n}^{x}-u_{l,m,n}^{x})-(u_{l+1,m,n}^{y}-u_{l,m,n}^{y})-(u_{l+1,m,n}^{z}-u_{l,m,n}^{z})\\ +(u_{l-1,m,n}^{x}-u_{l,m,n}^{x})-(u_{l-1,m,n}^{y}-u_{l,m,n}^{y})-(u_{l-1,m,n}^{z}-u_{l,m,n}^{z})\\ +(u_{l,m,n+1}^{x}-u_{l,m,n}^{x})+(u_{l,m,n+1}^{y}-u_{l,m,n}^{y})-(u_{l,m,n+1}^{z}-u_{l,m,n}^{z})\\ +(u_{l,m,n-1}^{x}-u_{l,m,n}^{x})+(u_{l,m,n-1}^{y}-u_{l,m,n}^{y})-(u_{l,m,n-1}^{z}-u_{l,m,n}^{z})\\ +(u_{l,m+1,n}^{x}-u_{l,m,n}^{x})-(u_{l,m+1,n}^{y}-u_{l,m,n}^{y})+(u_{l,m+1,n}^{z}-u_{l,m,n}^{z})\\ +(u_{l,m-1,n}^{x}-u_{l,m,n}^{x})-(u_{l,m-1,n}^{y}-u_{l,m,n}^{y})+(u_{l,m-1,n}^{z}-u_{l,m,n}^{z})\\ +(u_{l-1,m-1,n-1}^{x}-u_{l,m,n}^{x})+(u_{l-1,m-1,n-1}^{y}-u_{l,m,n}^{y})\\ +(u_{l-1,m-1,n-1}^{z}-u_{l,m,n}^{z})+(u_{l+1,m+1,n+1}^{x}-u_{l,m,n}^{x})\\ +(u_{l+1,m+1,n+1}^{y}-u_{l,m,n}^{y})+(u_{l+1,m+1,n+1}^{z}-u_{l,m,n}^{z})\Big]\\ +C_2(u_{l,m+1,n+1}^{x}-2u_{l,m,n}^{x}+u_{l,m-1,n-1}^{x}) \end{equation}

There are similar equations for the displacements in the $y-$ and $z-$directions. The normal modes must be eigenfunctions of the symmetry of the system. Because of the translational symmetry of the crystal, the normal modes must have the form,

\begin{equation} \vec{u}_{l,m,n}=\vec{u}_{\vec{k}}\exp\left( i(l\vec{k}\cdot\vec{a}_1+m\vec{k}\cdot\vec{a}_2+n\vec{k}\cdot\vec{a}_3-\omega t)\right). \end{equation}

Substituting this form for the normal modes into the equations of motion above results in a set of algebraic equations that can be written in matrix form,

$$\frac{M}{C_1}\omega^2\vec{u}_{\vec{k}}= \textbf{m} \vec{u}_{\vec{k}},$$

Where the elements of the matrix $\textbf{m}$ are:

$$m_{11} = \frac{2}{3}\left[-\cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) - \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right) +4\right]- 2\frac{C_2}{C_1}\left(\cos\left( ak_x)-1\right)\right)$$ $$m_{12} = m_{21} = \frac{2}{3}\left[\cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) + \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right)\right] $$ $$m_{13} = m_{31} = \frac{2}{3}\left[ \cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) +\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) - \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right) \right]$$ $$m_{22} = \frac{2}{3}\left[-\cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) - \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right) +4\right]- 2\frac{C_2}{C_1}\left(\cos\left( ak_y\right)-1\right)$$ $$m_{23} = m_{32} = \frac{2}{3}\left[ -\cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) +\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) + \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right) \right] $$ $$m_{33} = \frac{2}{3}\left[ -\cos\left( \frac{a}{2}(-k_x+k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y-k_z)\right) - \cos\left( \frac{a}{2}(k_x-k_y+k_z)\right) -\cos\left( \frac{a}{2}(k_x+k_y+k_z)\right) +4\right] - 2\frac{C_2}{C_1}\left(\cos\left( ak_z\right)-1\right)$$

The eigenvalues $\lambda = \frac{M\omega^2}{C_1}$ can be found by setting the determinant equal to zero, $|\textbf{m} -\lambda \textbf{I}| = 0$.

This can be written as a cubic equation in $\lambda$.

$$-\lambda^3 +(m_{11}+m_{22}+m_{33})\lambda^2+(m_{12}^2+m_{13}^2+m_{23}^2 - m_{11}m_{22} - m_{11}m_{33} - m_{22}m_{33})\lambda + m_{11}m_{22}m_{33} +2m_{12}m_{13}m_{23} - m_{12}^2m_{33} - m_{13}^2m_{22}-m_{23}^2m_{11} =0$$

Cubic equations can be solved using Cardano's formula. The standard form of a cubic equation is,

$$a\lambda^3 + b\lambda^2 + c\lambda + d =0$$

Code that will calculate the coefficients $a,b,c,d$ and the matrix elements from the $\vec{k}$-vector is,

$k_xa=$  $k_ya=$  $k_za=$   $\frac{C_2}{C_1} = $  

$\sqrt{\frac{M}{C_1}}\omega = $

$\Gamma: 0,0,0$

$H: 0,0,\frac{2\pi}{a}$

$P: \frac{\pi}{a},\frac{\pi}{a},\frac{\pi}{a}$

$N: 0,\frac{\pi}{a},\frac{\pi}{a}$