Question: Given a collection of point-particles, interacting through an attractive force $\sim \frac{1}{r^2}$.

Knowing only $m\_1a=\sum\_i \frac{Gm\_1m\_i}{r^2}$ and initial conditions we can deduce the motion of the system.

Consequently we can observe that three quantities remains constant
A) center of mass of the system
B) total energy
C) angular momentum

How can we derive these 3 facts directly from $m\_1a=\frac{Gm\_1m\_2}{r^2}$ ?

Are these quantities conserved for any attractive force $\sim\frac{1}{r^n}$ ?

Given any monotonically decreasing force for $r$ in $(0,\infty)$, which are the conserved quantities?

Answer: 
Actually, this is the basic stuff every mechanics textbook should have. 

1. Center of mass is the most basic it needs just Newton laws:  

Second: $m\_i\ddot{\vec{r}\_i}= \vec{F}\_i$  

And third: $\sum\_i\vec{F}\_i = 0$  

Summing over i one obtains:  

$\frac{d^2}{dt^2}\left(\sum\_i m\_ir\_i\right) = \sum\_i\vec{F}\_i = 0$
2. For the total energy you need those forces to be potential and independent on time:  

$F\_i = -\frac{d}{d\vec{r}\_i}U(\vec{r}\_1,\vec{r}\_2,...)$  

Then just take the total energy :
$E = \sum\_i \frac{m|\vec{r}\_i|^2}{2}+U$  

And differentiate with respect to time:  

$\frac{dE}{dt} = \sum\_i m\_i(\ddot{\vec{r}}\_i\dot{\vec{r}}\_i)+\sum\_{i} \frac{dU}{d\vec{ r}\_i}\dot{\vec{r}}\_i = 0$
3. Finally, for the angular momentum $U$ must be rotational invariant:  

$U(R\vec{r}\_1,R\vec{r}\_2,R\vec{r}\_3,...) = U(\vec{r}\_1,\vec{r}\_2,\vec{r}\_3,...)$  

where R is the rotation matrix. Consider now very small (infinitesimal) rotation:  

$R\vec{r}\_i = \vec{r}\_i+\left[\delta\vec\phi\times\dot{\vec{r}}\_i\right]$  

Substituting and expanding, one can get:  

$U+\delta\vec\phi\sum\_i\left[\frac{dU}{d\vec{r}\_i}\times\vec{r}\_i\right]=U$  

Which works for every angle $\delta\phi$, so the following must hold
$\sum\_i\left[\frac{dU}{d\vec{r}\_i}\times\vec{r}\_i\right]=0$ 

Finally take the angular momentum:  

$\vec{M} = \sum\_i m\_i \left[\vec{r}\_i\times\dot{\vec{r}}\_i\right]$  

And differentiate with respect to time:  

$\frac{d\vec{M}}{dt} = \sum\_i \left[\vec{r}\_i\times m\_i\ddot{\vec{r}}\_i\right]+
\sum\_i m\_i \left[\dot{\vec{r}}\_i\times \dot{\vec{r}}\_i\right] = \sum\_i\left[\frac{dU}{d\vec{r}\_i}\times\vec{r}\_i\right]=0$

Phew... You can see that I never used the "pairwise" interactions. And I only needed the basic assumptions about the potential energy. 

Concerning "extra" things that you cans say for your systems -- check [the Virial theorem](http://en.wikipedia.org/wiki/Virial_theorem)
All codes of the N-body simulation make use of instantaneous propagation of the speed of gravity.
The other answers follow the same criterion that is good enough for local simulations.

**As a consequence of the finite speed of propagation of gravity and the motion of the masses there is a dissipative component.**

As the bodies are never static, the masses continuously lose energy to space because there is no known mechanism that reverses the direction of flow.

(An equivalent effect applies to the radiation of charged particles by virtue of its continued accelerated motion.)

from **Jorgen Kalckar and Ole Ulfbeck experiment** (a) (energy of two masses and a spring):
...

> 
> The difference between these two
>  energies is lost by each mass: it is
>  taken away by space-time, in other
>  words, it is radiated away as
>  gravitational radiation
> 
> 
> gravity can behave like a wave:
>  gravity can radiate
>  ...
>  All this follows from the expression of universal gravity when applied to moving observers, with the requirement that neither observers
>  nor energy can move faster than c.
> 
> 
> 

(a) ***motionmountain*** free online book, ch 18, pag 526/527: Motion in General Relativity, Gravitational Waves.