In the absence of external forces, the motion of the centre of mass of the whole system will be a straight line. Total momentum must be conserved, and this may be written as the total mass times the centre-of-mass velocity vector; or, equivalently, as the sum of the centre-of-mass momenta of the two composite rigid bodies:
$$
M \dot{\mathbf{R}} =
\sum_j m_j \dot{\mathbf{r}_j} =
\sum_{j\in 1} m_j \dot{\mathbf{r}_j} + \sum_{j\in 2} m_j \dot{\mathbf{r}_j} =
M_1\dot{\mathbf{R}}_1 + M_2\dot{\mathbf{R}}_2
= \text{constant}
$$
where $M=\sum_j m_j$, $M_1=\sum_{j\in 1} m_j$, $M_2=\sum_{j\in 2} m_j$, and we have defined the centre of mass positions
$$
\mathbf{R} = \frac{1}{M}\sum_j m_j \mathbf{r}_j , \quad
\mathbf{R}_1 = \frac{1}{M_1}\sum_{j\in 1} m_j \mathbf{r}_j , \quad
\mathbf{R}_2 = \frac{1}{M_2}\sum_{j\in 2} m_j \mathbf{r}_j .
$$
Most likely the problem in your simulation arises from the way you have implemented the assumption you made in your last paragraph, since the gravitational interaction actually acts between each pair of point particles. Having summed up all these forces acting on body 1, this determines the acceleration of the centre of mass of body 1; similarly for body 2. But (by Newton's third law) there will be no effect on the total (centre of mass) momentum of the whole system. Also, you would need to work out the total torque acting on body 1 and consider its rotation about the centre of mass; similarly for body 2.
(If you have neglected to consider torques, and the rotation of the two rigid bodies, this is somewhat unphysical, but has no bearing on the underlying issue of the motion of the two centres of mass).
