Your force is correct, that is also the expression for $\ddot{r}$ in the real Schwarzschild De Sitter metric when you set the first proper time derivatives of the spatial coordinates to zero:
The geodesic equation gives the radial component of the 4-acceleration (in natural units):
$$ \ddot{r}= \color{gray}{ \frac{\left(\Lambda r^3-3\right) \dot{r}^2}{r \left(\Lambda r^3-3 r+6\right)} } -\frac{\left(\Lambda r^3-3 r+6\right) \left(\color{gray}{ 3 r^3 \left(\dot{\theta}^2 +\sin ^2 \theta \ \dot{\phi}^2\right) }+\left(\Lambda r^3-3\right) \dot{t}^2\right)}{9 r^3} $$
where you set $\dot{r}=\dot{\rm \theta}=\dot{\rm \phi}=0$ and plug in
$$ \dot{t}=\sqrt{g^{t t}} \ \color{gray}{ \gamma } = \sqrt{\frac{1 \ / \ (1-2/r-\Lambda r^2/3)}{1-\color{gray}{ v^2 }}}$$
with $v=0$, where $v$ is the velocity measured by local and stationary (relative to the dominant mass) Fidos, then you get
$$\ddot{r} = -\frac{1}{r^2}+\frac{\Lambda r}{3}$$
which is, in natural units, the expression you correctly guessed. The overdot is the differentiation with respect to proper time, but in the newtonian limit the proper and coordinate time are the same.
This equation assumes the dominant mass in the center, for an n-body simulation the $r$ in the $-1/r^2$ term is relative to the masses and the $r$ in the $+\Lambda r/3$ term relative to the center of your coordinate grid from which everything else accelerates away.
