I am considering this more of a math question than a physics question. Note that gravity can be replaced by other interactions that follow an inverse-square law in 3D Euclidean space. BTW, here is a similar question, but not quite similar enough that I should stay quiet, IMO: Can I assign a gravity field to an infinite grid of point masses?
Consider an $s*s*s$ cube. Each side is like a portal to the opposite side, conserving rotation and direction.
One way to view this space is to copy-paste the cube to an infinite 3D grid surrounding it, thus obtaining an Euclidean space. OK, but if there are a number of particles in the wraparound space, the corresponding Euclidean space will have an infinity of particles, possibly making it less computationally tractable.
Anyhow, I have found the Euclidean conversion to be a useful instrument for trying to understand this. Take the case where there are two particles $p_1,p_2$. (The general n particle case can be computed by superposition of scaled and shifted simple cases). $p_1$ is in the center (origo) of the cube, and $p_2$ in position $\vec r$, where $|r_x|,|r_y|,|r_x| \le s/2$. The objective is to calculate the net gravity working on $p_1$. One observation I make, is that for every copy of $p_1$ there is an opposing copy particle which pulls with the same force in the opposite direction. Therefore $p_1$ does not pull on itself. (Likewise, $p_2$ does not pull itself.)
The contribution from $p_2$ is an infinite sum over the grid. Something like (hope pseudocode is ok):
sum := 0 //vector
for i in the set of integers:
for j in the set of integers:
for k in the set of integers:
rijk = s*(i,j,k) + r //vector
sum := sum + C*(rijk/|rijk|)/|rijk|^2
//Or equivalently: sum := sum + C*(rijk/|rijk|^3)
Output sum //vector holding the force acting by p2 on p1
The relative contribution of the base r term gets smaller the bigger i,j,k get, so I expect that for bigger i,j,k the opposing -i,-j,-k cell will have a very similar impact in the opposite direction, but sadly the number of cells around distance $r_{ijk}$ scales with $r_{ijk} ^2$, so AFAICS even small differences can add up to something significant. If it can be proved that the sum eventually converges for increasing $r_{ijk}$, the computation can be approximated by a finite set of cells, and this will actually be a quite satisfying solution. But a proof is needed.
But maybe there is a better way to compute this than using the transformation to Euclidean space. Maybe applying the inverse square law is not necessary at all, because that is something which belongs in 3D Euclidean space. Maybe the considerations that give rise to the inverse square laws in Euclidean space can be put better to use in the 3-torus space? I am very open for suggestions. :-)
An application for knowledge on this type of interaction is to do approximate simulations of small physical volumes of homogenous material (ignoring quantum effects for now). I can model electric and gravitational interactions between a few hundred or a few thousand particles in an open space or a closed box (with wall bouncing), but they are poor approximations for a small volume in a homogenous material. I am thinking that a wrapped space gives a better approximation of the microvolume's interaction with its environment.