In the wikipedia page on BPS bound (https://en.wikipedia.org/wiki/Bogomol%27nyi%E2%80%93Prasad%E2%80%93Sommerfield_bound) there is a proof, but I do not understand the final step:
$$E\geq \pm \dfrac 1g\int d^3x \text{Tr}\left[\vec{D\varphi}\cdot \vec B\right]=\pm \dfrac 1g\int_{S^2\text{boundary}} \text{Tr}\left[\varphi\vec B\cdot d\vec S\right]$$
Here $D$ is the covariant derivative. I guessed there is an integration by parts, so there is
$$\int d^3x \text{Tr}\left[\vec{D\varphi}\cdot \vec B\right]= \int d^3x \text{Tr}\left[(D_i\varphi) B_i\right] = \int d^3x \text{Tr}\left[D_i(\varphi B_i)\right] - \int d^3x \text{Tr}\left[\varphi D_i( B_i)\right]$$ By divergence theorem the first term becomes $$\int d^3x \text{Tr}\left[D_i(\varphi B_i)\right] = \int_{S^2\text{boundary}} \text{Tr}\left[\varphi\vec B\cdot d\vec S\right]$$
But what about the second term? Why does it disappear?
In this case there is an assumption that the magnetic monopole exists, so definitely this is not applying $\vec D \cdot \vec B=0$.
In Yang-Mills theory, the Jacobi identity, reflecting associativity of covariant derivatives, is but the full antisymmetrization of three covariant derivatives, so, by 7.24, or, in components, we have $$ D_\mu \tilde G ^{\mu \nu}=\partial_\mu \tilde G ^{\mu \nu} +[A_\mu, \tilde G ^{\mu \nu}] \\ = \epsilon^{\mu\nu\rho\kappa} \Bigl ( \partial_\mu (\partial_\rho A_\kappa + A_\rho A_\kappa )+[A_\mu, (\partial_\rho A_\kappa + A_\rho A_\kappa )] \Bigr ) =0, $$ where, of course, $$\tilde G ^{\mu \nu} \equiv \epsilon^{\mu\nu\rho\kappa} G _{\rho \kappa}/2$$ is the dual field strength and we have used the antisymmetric properties of the totally antisymmetric tensor $\epsilon^{\mu\nu\rho\kappa}$. This is the fundamental YM Bianchi identity.
When the index $\nu$ is 0, all other indices are spacelike, and the dual field strength reduces to a magnetic field.
Proceed to check that, in sharp contrast to your confused chain rule misimpression, $$ \operatorname{Tr} \partial_k (\epsilon_{ijk} G_{ij} \phi)=\operatorname{Tr} \epsilon_{ijk} (D_k G_{ij} \phi+ G_{ij} D_k \phi) . $$
Do you now see how the covariant completion terms in the covariant derivatives cancel each other, given the cyclicity of the trace ? The total divergence had better not be a gauge-covariant divergence!
The first term on the r.h.s. then, proportional to $D_k \tilde G ^{k0}$ vanishes by the above identity, which leaves you with just the second term, as in the WP sound proof.