I found the following equations in a dynamics textbook,

The gravitational potential energy for any two particles in a n-partice system is given by,

$$V_{ij} = - \frac {G m_i m_j}{r_{ij}}$$

where $$r_{ij}$$ is the distance between $$m_i$$ and $$m_j$$. The total potential energy of the system is,

$$V = \frac{1}{2} \sum_{i = 1}^{n} \sum_{j = 1}^{n}V_{ij}$$ (i $$\neq$$ j)

If $$R_i$$ is the postion vector of the $$i^{th}$$ particle, then

$$\frac{\partial{V}} {\partial {\vec{R_{i}}}} = - \frac{\partial{V}}{\partial{\vec{r_{ji}}}} + \frac{\partial{V}}{\partial{\vec{r_{ij}}}} = -2 \frac{\partial{V}}{\partial{\vec{r_{ij}}}}$$

What does it mean to take the derivative of a Scalar Function$$(V)$$ with respect to a vector$$(\vec{R_1})$$? Is it directional derivative?

I've been trying all day to get the last equation. I would be very grateful if somebody could help me(or mention some reference perhaps). I don't really know which part of math is used to get the last equation.

Let’s do it for two particles – generalization will be straightforward. The potential energy of two particles is $$V_{ij} = - \frac {Gm_i m_j}{r_{ij}}=-\frac {Gm_i m_j}{|\vec{r_{ij}}|}=- \frac {Gm_i m_j}{|\vec{R_{i}}-\vec{R_{j}}|}$$, because $$r_{ij} =|\vec{R_{i}}-\vec{R_{j}}|$$.
By definition $$\frac{\partial}{\partial\vec{R}}V({R_{i}})=\frac{\partial} {\partial{R_x}}\vec{e_{x}}+\frac{\partial}{\partial{R_y}}\vec{e_y}+\frac{\partial} {\partial{R_z}}\vec{e_z}$$, where vectors $$\vec{e_x}, \vec{e_y}, \vec{e_z}$$ are directed along axis X, Y and Z correspondingly and each has a unit length.
Next, $$\frac{\partial}{\partial\vec{R_1}}\frac{1}{|\vec{R_1}-\vec{R_2}|}=(\frac{\partial} {\partial{R_{1x}}}\vec{e_x}+\frac{\partial}{\partial{R_{1y}}}\vec{e_y}+\frac{\partial} {\partial{R_{1z}}}\vec{e_z})\frac{1}{\sqrt{(R_{1x}- R_{2x})^2+(R_{1y}- R_{2y})^2+(R_{1z}- R_{2z})^2}}$$.
Taking derivatives we get explicitly: $$\frac{\partial}{\partial\vec{R_1}}\frac{1}{|\vec{R_1}-\vec{R_2}|}=-\left((R_{1x}- R_{2x})\vec{e_x}+(R_{1y}-R_{2y})\vec{e_y}+(R_{1z}- R_{2z})\vec{e_z}\right)(\frac{1}{\sqrt{(R_{1x}- R_{2x})^2+(R_{1y}- R_{2y})^2+(R_{1z}- R_{2z})^2)}})^{3}=-\vec{r_{12}}\left((R_{1x}- R_{2x})^2+(R_{1y}- R_{2y})^2+(R_{1z}- R_{2z})^2\right)^{-\frac{3}{2}}$$
If you take the derivatives with respect to $$\vec{R_2}$$ you will get the same formula, but with the negative sing.