Note: I am working in the Lorentz-Heaviside system and all the integrals are over the whole space.
Definitions: $$\vec E= \vec E_f+\vec E_b$$ $$\phi=\phi_f+\phi_b$$ $$\vec D=\vec E+\vec P$$ $$\rho=\rho_f+\rho_b$$ Here the subscripts f and b refer to free and bound charges, thus $\vec E_f$ is the field created by free charge when all the charges are at the desired (final) positions.
$\\$
The total energy of the system is $$W=\frac12\int\rho\phi\, dV=\frac12\int E^2\, dV$$
But we are only interested in the energy that is required to move free charge from infinity to some desired position, thus we need to subtract the energy needed to create the system of the bound charge itself, namely $$\frac12\int \rho_b\phi_b\,dV$$ Since this is the energy of creation of bound charge system, where $\phi_b$ is the final potential caused by this system. One may think to also subtract the term of $\rho_b\phi_f\equiv\rho_f\phi_b$ but this is not correct since this is the interaction energy of two systems and must be included. Thus we get $$W=\frac12\int\rho\phi\, dV-\frac12\int \rho_b\phi_b\,dV$$ $$W=\frac12\int(\rho_f+\rho_b)(\phi_f+\phi_b)\, dV -\frac12\int \rho_b\phi_b\,dV$$ $$W=\frac12\int\rho_f\phi_f\, dV+\int\rho_f\phi_b\, dV\tag{1}$$ $$W=\frac12\int\rho_f(\phi+\phi_b)\, dV$$ $$W=\frac12\int\vec E\cdot\vec D\, dV+\frac12\int\rho_f\phi_b\, dV\tag{2}$$
$\\$
We can also derive this as follows: Consider that initially when the free charge is at infinity, we bring $dq_f$, the work needed for this will be $$dw_1=dq_f\int_{\phi_{b_1}}^{\phi_{b_2}}\phi_b$$ For the second free charge $$dw_2=dq_f\int_{\phi_{b_2}}^{\phi_{b_3}}\phi_b+dq_f\int_{\phi_{f_1}}^{\phi_{f_2}}\phi_f$$ Where the second integral comes due to the first charge.
Following the pattern, we get the net work to be $$W=\int\rho_f\phi_b\, dV+\frac12\int\rho_f\phi_f\, dV$$ Here both the potentials are the final potentials of the respective configurations. This is the same as equation 1.
Thus the question is why are we getting an extra integral in equation 2? What is going wrong in this derivation?