How to calculate magnetic force using differentiation of magnetic energy?

I tried to calculate the magnetic force between two circular current loops using numeric integration and differentiation of magnetic energy: $$\int d^3\vec{r} \frac{1}{2 \mu_0} \vec{B}^2 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(1).$$ But after repeated examination, the result still comes out to be essentially exactly the negative of what I calculated using the force law $$\vec{F}= \int I d\vec{l} \times \vec{B} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(2) ,$$ where the magnetic field is got using the Biot–Savart law $$\vec{B}(\vec{r})=\frac{\mu_0}{4 \pi } \int d^3 \vec{r}' \frac{I d\vec{l} \times \vec{r}' }{|\vec{r}'|^3} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(3) .$$

I mean, when I calculate the total magnetic energy of two circular loops with opposite current in a configuration like in the following figure, they should repel right? And thus the longer the distance between them the lower the total magnetic energy right? But my result is the longer the distance between them the higher the total magnetic energy. But what is even more strange is that the differentiation of the numerically calculated magnetic energy essentially exactly matches in magnitude the force calculated using the force law (which, of course, give a repulsive result).

I wonder why, or would someone please redo this calculation and see if I get it right or wrong.

For convenience of possible recalculation, the following analytical formula may be useful.

For a circular loop with current I, radius R, the result is, $$B_z(\rho, z)= \frac{\mu_0 I}{2 \pi} \frac{1}{\sqrt{(\rho+R)^2+z^2}} [K(k'^2)+ \frac{R^2-\rho^2-z^2}{(R-\rho)^2+z^2} E(k'^2)],\,\,\,\,\,\,\,\,(4)$$. $$B_r(\rho, z)= \frac{\mu_0 I}{2 \pi} \frac{z}{\rho\sqrt{(\rho+R)^2+z^2}} [-K(k'^2)+ \frac{R^2+\rho^2+z^2}{(R-\rho)^2+z^2} E(k'^2)] ,\,\,\,\,\,\,\,\,(5)$$ $$k'=\sqrt{\frac{4 R \rho}{(\rho+R)^2+z^2}} ,\,\,\,\,\,\,\,\,(6)$$ I got this from https://summit.sfu.ca/_flysystem/fedora/2023-04/DanielCarletonFinalUGThesisEngSci2022.pdf except for some minor corrections. There seems to be some typo there. For equation (4) - (6) which I use I confirmed using integration using the Biot–Savart law equation (3).

This expression has singularity at the wire position ($$\rho=R, z=0$$). When using this to calculate magnetic energy (1) by integrations, the singularity remains. I made some handwavy estimation that when calculating the differentiation of the magnetic energy to get the force such singular regions actually do not give finite contribution. So when I actually perform the calculations, I simply exclude spartial regions within distance $$\epsilon$$ from the wire. $$\epsilon$$ is chosen to be much smaller than other length scales in the question. Since the final result matches the result got by the force law in magnitude so well, I tend to guess this estimation is valid.