I have angular momenta $S=\frac{1}{2}$ for spin, and $I=\frac{1}{2}$ for nuclear angular momentum, which I want to add using the Clebsch-Gordan basis, so the conversion looks like:
$$ \begin{align} \lvert 1,1\rangle &= \lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr)\tfrac{1}{2}\tfrac{1}{2}\rangle,\tag{4.21a} \\ \lvert 1,0\rangle &= \frac{1}{\sqrt{2}}\biggl(\lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr)\tfrac{1}{2},-\tfrac{1}{2}\rangle + \lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr),-\tfrac{1}{2}\tfrac{1}{2}\rangle\biggr),\tag{4.21b} \\ \lvert 1,-1\rangle &= \lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr),-\tfrac{1}{2},-\tfrac{1}{2}\rangle,\tag{4.21c} \\ \lvert 0,0\rangle &= \frac{1}{\sqrt{2}}\biggl(\lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr)\tfrac{1}{2},-\tfrac{1}{2}\rangle - \lvert\bigl(\tfrac{1}{2}\tfrac{1}{2}\bigr),-\tfrac{1}{2}\tfrac{1}{2}\rangle\biggr),\tag{4.21d} \end{align} $$
where $F=I+S$, so this is the basis $\lvert F m_F \rangle = \sum_m \lvert\bigl(I S\bigr),m_I m_S\rangle $.
Now since adding angular momenta is commutative, the exchange between $I$ and $S$ shouldn't mathematically introduce any kind of difference.
In other words, in the basis described in those equations, it shouldn't matter whether I write it as $\lvert\bigl(I S\bigr),m_I m_S\rangle$ or $\lvert\bigl(S I\bigr),m_S m_I\rangle$, right?
Now the problem is the following: I have created the Hamiltonian matrix $H=-\vec{\mu}\cdot \vec{B} = -2 \mu B_z S_z/\hbar$ in the $\lvert F m_F \rangle$ representation, and actually the result depends on how you call those angular momenta, so the result could be
$$H = \begin{pmatrix} \mu_B B & 0 & 0 & 0 \\ 0 & - \mu_B B & 0 & 0 \\ 0 & 0 & 0 &\mu_B B \\ 0 & 0 & \mu_B B & 0 \end{pmatrix}$$
Or could be
$$H = \begin{pmatrix} \mu_B B & 0 & 0 & 0 \\ 0 & - \mu_B B & 0 & 0 \\ 0 & 0 & 0 &-\mu_B B \\ 0 & 0 & -\mu_B B & 0 \end{pmatrix}$$
Depending on how you "label" them, $I$ or $S$... which is very confusing!
This happens because the off-diagonal terms
$$\left\langle 1 0 \right| S_z \left| 0 0 \right\rangle = \frac{1}{2} \left( \left\langle (\frac{1}{2} \frac{1}{2}) \frac{1}{2} -\frac{1}{2} \right| + \left\langle (\frac{1}{2} \frac{1}{2}) -\frac{1}{2} \frac{1}{2} \right| \right) S_z \left( \left| (\frac{1}{2} \frac{1}{2}) \frac{1}{2} -\frac{1}{2} \right\rangle - \left| (\frac{1}{2} \frac{1}{2}) -\frac{1}{2} \frac{1}{2} \right\rangle \right)$$
will be either $\hbar/2$ or $-\hbar/2$ depending on your convention whether it's $\lvert\bigl(I S\bigr),m_I m_S\rangle$ or $\lvert\bigl(S I\bigr),m_S m_I\rangle$.
How can I understand this physically and mathematically? Shouldn't the addition be commutative and the process be blind to which labels I use?