It is well known that, if I have a system of $N$ particles acted upon only by conservative internal and external forces, then I can obtain the force on the $\mathrm{i^{th}}$ particle as $$\textbf{F}_i = -\nabla_iU$$ where the potential energy of the entire system $U= U_{int} + U_{ext}$ is the sum of internal (due to inter-particle interactions) and external (due to externally impressed forces) potential energies. Now suppose I want the net force $\textbf{F}$ on the system of particles. My first question is: I suppose this is $\textit{defined}$ as the sum of forces from external sources on the individual particles? If so, then I should be able to proceed, assuming that the internal forces obey Newton's third law, with $$\sum_{i}-\nabla_i U=\sum_{i}-\nabla_i U_{ext} = \sum_{i}\textbf{F}_{ext} \equiv \textbf{F},$$ where in the first equality I've observed that internal forces cancel pairwise given the obeyance of Newton's third law. Now this entire setup has been for the following question, where my goal is to obtain the Coulomb force as a theorem from an assumed form for the Coulomb potential energy.
Suppose that a finite charge distribution $\rho$ is of the form $\rho(\textbf{r}) = \rho_2(\textbf{r}-\textbf{R})$ where $\textbf{R}$ is a parameter "locating" the body and, per this answer, would usefully be placed at the centre of mass of the distribution if we were interested in the dynamics of the system. Plunging forward, we ask about the force on $\rho$ due to an external electrostatic potential $\varphi_1$, where we suppose that the potential energy $U$ of the system is given by (an integral with no bounds is assumed to be "over all space") $$U \equiv \int d^3r \ \rho(\textbf{r})\varphi_1(\textbf{r})=\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}).$$ Now I want to prove the following theorem for this particular electrostatic system: $\textbf{F} = -\nabla_\textbf{R}U$ where, importantly, $\textbf{F}$ takes the typical Coulomb form. Is the following correct? $$-\nabla_\textbf{R}U = -\nabla_\textbf{R}\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r})=\int d^3r \ (-\nabla_\textbf{R})\rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}) \stackrel{(1)}{=} \int d^3r \ \nabla_\textbf{r}\rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}) \stackrel{(2)}{=} \\-\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\nabla_\textbf{r}\varphi_1(\textbf{r} )\stackrel{(3)}{=} \int d^3r \ \rho_2(\textbf{r}-\textbf{R})\textbf{E}_1(\textbf{r}) = \int d^3r \ \rho(\textbf{r})\textbf{E}_1(\textbf{r}) \stackrel{(4)}{=} \textbf{F}$$ where in (1) I've changed the variables with respect to which the gradient is taken, in (2) I've used an integration by parts and that the integrand vanishes on a boundary at infinity (since the charge distribution is finite), in (3) I've used the definition of the electric field as the gradient of its potential, and in (4) I've used the identification of $\textbf{F}$ as the sum of external forces (in the continuum limit) as I asked earlier in the question. Now it seems like what I've done intimately requires that the charge distribution be finite (in particular, in the aforementioned integration by parts step), whereas the second equation in this post seems to suggest that my result should hold in the framework of classical mechanics no matter the nature of $\rho$. Have I made a mistake somewhere then in either analysis?