In statistical physics, one finds the probability distribution $\rho[q]$ that maximizes the entropy $S$
$$ S=-k_B \sum_{q\in\mathbb{Q}} \rho[q]\ln \rho[q]\tag{1} $$
under the constraints:
$$ \begin{align} \overline{E}&=\sum_{q\in\mathbb{Q}}E[q]\rho[q] \tag{2}\\ 1&=\sum_{q\in\mathbb{Q}}\rho[q] \tag{3} \end{align} $$
by using the method of the Lagrange multipliers:
$$ \mathcal{L}= \left( \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left(-1+ \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( - \overline{E} + \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) \tag{4} $$
To obtain the Gibbs measure, one must find the stationary action
$$ \frac{\partial \mathcal{L}}{\partial \rho[q]}=0 $$
I would like to know if my approach to it is valid. Let me start:
$$ \begin{align} &\frac{\partial }{\partial \rho[q]} \left[ \left( -k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left(-1+ \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( - \overline{E} + \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) \right] =0 \\ &\frac{\partial }{\partial \rho[q]} \left[ \left( -k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left( \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) \right] =0 \\ & \frac{\partial }{\partial \rho[q]} \left( k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) = \frac{\partial }{\partial \rho[q]} \left[ \lambda \left( \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) \right]\\ & \left( k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) = \lambda \left( \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) +C\\ & S = \lambda +C + \beta \overline{E} \\ & S = \ln Z + \beta \overline{E} \end{align} $$
where we define $\ln Z = \lambda+C$.
I struggle with my derivation because it seems like there is an arbitrary choice to make. For instance, if I would have instead done this:
$$ \begin{align} &\frac{\partial }{\partial \rho[q]} \left[ \left( -k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left(-1+ \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( - \overline{E} + \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) \right] =0 \\ &\left( -k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left(-1+ \sum_{q\in\mathbb{Q}} p[q] \right) + \beta \left( - \overline{E} + \sum_{q \in \mathbb{Q}} E[q] \rho[q] \right) +C =0 \\ &\left( -k_B \sum_{q\in\mathbb{Q}} \rho[q] \ln \rho[q] \right) + \lambda \left(-1+ 1 \right) + \beta \left( - \overline{E} + \overline{E} \right) +C =0 \\ &S + C =0 \end{align} $$
Why different results if I simply change the order in which I make the replacements? Which, if any, is correct?