Unitary matrix that transforms Bloch wavefunctions of k-points connected by crystal symmetries In the following link https://docs.abinit.org/theory/wavefunctions/#symmetry-properties, there's Eqn. (8) that describes for degenerate states how wavefunctions of k-points connected by symmetry $\mathcal R$ are connected (transformed from each other):
$$ \begin{equation} \Psi_{n\bf k} \left( \mathcal R^{-1} \left( \bf r-t \right)  \right) = \sum_{\alpha \in \mathcal C_{n\bf k}} D_{\alpha n}\left( \mathcal R \right) \Psi_{\alpha\mathcal R \bf k} \left( \bf r \right) \end{equation} $$
My question is, how does one calculate or obtain the matrix elements $D_{\alpha n}$?
 A: Let's just take your equation seriously. We have
$$ 
\Psi_{n\bf k} \left( \mathcal R^{-1} \left( \bf r-t \right)  \right) = \sum_{\alpha \in \mathcal C_{n\bf k}} D_{\alpha n}\left( \mathcal R \right) \Psi_{\alpha\mathcal R \bf k} \left( \bf r \right)
$$
This tells us that $\Psi_{n\bf k} \left( \mathcal R^{-1} \left( \bf r-t \right)  \right)$ can be written as a linear combination of $\Psi_{\alpha\mathcal R \bf k} \left( \bf r \right)$. If we want to know the coefficients of the linear combination, we should just multiply both sides of your equation by $\Psi^*_{\beta\mathcal R \bf k} \left( \bf r \right)$ and integrate with respect to $\mathbf{r}$.
\begin{align*}
\int d\mathbf{r}\ \Psi^*_{\beta \mathcal{R}\mathbf{k}}(\mathbf r)\Psi_{n\bf k} \left( \mathcal R^{-1} \left( \bf r-t \right)  \right) &= \sum_{\alpha \in \mathcal C_{n\bf k}} D_{\alpha n}\left( \mathcal R \right) \int d\mathbf{r}\ \Psi^*_{\beta\mathcal R \bf k} \left( \bf r \right)\Psi_{\alpha\mathcal R \bf k} \left( \bf r \right)\\
&= \sum_{\alpha \in \mathcal C_{n\bf k}} D_{\alpha n}\left( \mathcal R \right) \delta_{\alpha\beta} \qquad\text{(because of orthogonality)}\\
&= D_{\beta n}\left( \mathcal R \right) \qquad\text{(using the $\delta$ to do the sum)} 
\end{align*}
So we have our formula for $D_{\beta n}$, as desired!
