What does it mean to take the “ gradient with respect to the position $r_ i$”? Let’s say we have a number of particles (charged, massive or anything that can create potential energy). The total potential energy of any particle can be given by $$U_i (\vec r_1, \vec r_2, ... \vec r_N) =  \sum_{j\gt i}^{N} U_{ij}(|r_i -r_j|)$$ and hence total potential of the system is 
$$
U_{total}(\vec r_1, \vec r_2, ..., \vec r_N)= \sum_{i=1}^{N} U_i$$
($U_{ij}$ means potential energy of the particle $i$ due to the particle $j$)
Now, to get the force on the particle $i$ we write 
$$
\mathbf F_i = -\nabla_{r_i} U_{total}(\vec r_1, \vec r_2, ..., \vec r_n)
$$
Now, that subscript after $\nabla$ is posing some problems for me. 


*

*What is meant by “gradient with respect to $r_i$? Isn’t $r_i$ a fixed position of $i$ th particle? It’s same as saying “take the derivative of $f$ with respect to number 2”. 

*Why not to take the simple gradient, that is only $\nabla$, of $U_{total} (\vec r_1, \vec r_2, ... \vec r_N)$ and evaluate the gradient at $r_i$ ? 
Can you please resolve my doubts? 
Here is a Stackexchange post which address the similar problem, but I’m quite unsatisfied with the answers over there. 
 A: The total potential energy is a function $U: \left(\Bbb{R}^3\right)^N \to \Bbb{R}$. Where the physical idea is that the $i^{th}$ copy of $\Bbb{R}^3$ tells us the position $(x_i, y_i, z_i)$ of the $i^{th}$ particle. So, writing something like $\mathbf{F}_i = - \nabla_{\mathbf{r}_i}U$ means: $\mathbf{F}_i : \left(\Bbb{R}^3\right)^N \to \Bbb{R}^3$ is the vector valued function defined as
\begin{align}
\mathbf{F}_i &= - \left( \dfrac{\partial U}{\partial x_i}\, \mathbf{e}_1 + \dfrac{\partial U}{\partial y_i}\, \mathbf{e}_2 + \dfrac{\partial U}{\partial z_i}\, \mathbf{e}_3\right) \\
&\equiv - \left(\dfrac{\partial U}{\partial x_i}, \dfrac{\partial U}{\partial y_i}, \dfrac{\partial U}{\partial z_i} \right),
\end{align}
where I use $\equiv$ to mean they're the same thing, expressed in different notation.
In other words, the force $\mathbf{F}_i$ on the $i^{th}$ particle is obtained by differentiating the total potential energy with respect to the $3$ cartesian coordinates of that particle.

Edit:
The short answer to the question in the comments is that "no, $\mathbf{F}_a = - \nabla_{\mathbf{r}_a}(U_{\text{total}})$" is a mathematical statement which needs to be proven and is not a "law". Here, we crucially make use of the fact that the potentials $U_{ij}$ depend only on $|\mathbf{r}_i - \mathbf{r}_j|$.
By definition, the total force on $a^{th}$ particle, due to all other particles is
\begin{align}
\mathbf{F}_a = \sum_{j \neq a} \mathbf{F}_{\text{$j$ on $a$}}
\end{align}
What is $\mathbf{F}_{\text{$j$ on $a$}}$? Well, you simply take the potential energy $U_{aj}$ and differentiate with respect to $\mathbf{r}_a$, i.e $\mathbf{F}_{\text{$j$ on $a$ }} = - \nabla_{a}(U_{aj})$ (for ease of typing, I use $\nabla_a$ to mean $\nabla_{\mathbf{r_a}}$, which I defined above to mean differentiation with respect to $x_a, y_a, z_a$). Hence,
\begin{align}
\mathbf{F}_a = \sum_{j \neq a} \mathbf{F}_{\text{$j$ on $a$}} = - \sum_{j \neq a} \nabla_{a}(U_{aj}) \tag{$*$}
\end{align}
I now claim that this is also equal to $-\nabla_a U_{\text{total}}$. To prove this, note first of all that $U_{ij} = U_{ji}$ and that $\nabla_{a}(U_{ij}) = 0$ if $a \notin \{i,j\}$ (because the potential energy depends only on the distance between the $i$ and $j$ particles, so clearly it doesn't depend on some other $\mathbf{r}_a$).
So, we have
\begin{align}
\nabla_a(U_{\text{total}}) &= \nabla_a\left( \sum_{i=1}^n \sum_{j>i} U_{ij}\right)
\end{align}
Now, we shall split the sum over $i$ into 3 pieces: $i=a$, $i<a$ and $i>a$. Then, we get:
\begin{align}
\nabla_a(U_{\text{total}}) &= \sum_{j>a} \nabla_a(U_{aj}) + \sum_{i<a} \sum_{j>i} \nabla_a(U_{ij}) +  \sum_{i>a} \sum_{j>i} \nabla_a(U_{ij})
\end{align}
Now, recall the above mentioned property that $\nabla_a(U_{ij}) = 0$ if $a \notin \{i,j\}$. So, the only non-zero terms in the above summation are:
\begin{align}
\nabla_a(U_{\text{total}}) &= \sum_{j>a} \nabla_a(U_{aj}) + \sum_{i<a} \nabla_a(U_{ia}) +  0 \\
&= \sum_{j \neq a} \nabla_a(U_{aj}), \tag{$**$}
\end{align}
where in the last line I used the fact that $U_{ia} = U_{ai}$, and I renamed summation indices and combined everything. So, if you combine $(*)$ and $(**)$ then you immediately see that
\begin{align}
\mathbf{F}_a &= - \nabla_a(U_{\text{total}}).
\end{align}
A: What is meant by taking the gradient with respect to $\vec{r}_i$ is to calculate the quantity
$$\vec{\gamma} = \left(\frac{\partial U_{total}}{\partial r_{i, 1}}, \frac{\partial U_{total}}{\partial r_{i, 2}}, \frac{\partial U_{total}}{\partial r_{i, 3}}\right)$$
where $r_{i, j}$ is the $j$th coordinate of the vector $\vec{r}_i$. What is crucial to note here is that neither the vector $\vec{r}_i$ nor its components are fixed quantities. They are variables, or degrees of freedom, that the system depends upon. That is why taking the derivative of a function with respect to them makes sense in the first place. 
