Let me suggest an alternative approach, whereby you model more from first principles, and end up not even needing to consider torque, etc, per se. They just automatically emerge from the first-principles model dynamics, as follows.
I'd imagine you're thinking of a dipole as two point charges $+q,-q$ subject to the constraint that they're separated by a constant distance, as though connected by a massless "virtual rigid rod" of length $R$ (as labelled in your diagram). And now you've got to think about torques, centers of mass, etc.
Now, since this is just a project (i.e., no buildings are going to fall down, etc), think about dipoles the following way, instead. You've still got your $+q,-q$ charges, as before, but they're no longer constrained by a "virtual rod", per se. Instead, they're subject to an additional rod-like-central-force of the form $F_{rod}=-k(d_q-R)^3$ that acts along the line joining $+q,-q$, where $d_q$ is the distance between $+q$ and $-q$, with $F$ attractive when $d_q>R$ (and repulsive when $d_q<R$). And $k$ is arbitrarily chosen very large, so that the "rod" is very "stiff", i.e., when $d_q\not=R$ the restoring force drags the two particles back to equilibrium $R$ pretty much immediately. And, of course, that $(d_q-R)^3$ is also arbitrary -- you could use some exponential function of $d_q$ and $R$, or any kind of function you like.
So now your entire programming model is just a bunch of central forces acting on point masses. You've presumably already programmed your Coulomb forces. Now, in addition, just program your $F_{rod}$ forces. Moreover, you described your model as an "ensemble of dipoles". Forget that description, per se. Consider your model described as an "ensemble of point charges" plus an "ensemble of edges", giving you a disconnected graph. Each "edge" connects the "vertices"/charges comprising a dipole. And the "edge weight" (are you familiar with all this from graph theory?) is your $R$, which can be different for each dipole (you could even parameterize "edges" to have different $k$'s for each dipole, and/or different $F_{rod}$ force laws, describing different families of charged particles).
So now you've got simple data structures with simple loops to sum up all the central forces (both Coulomb and our $F_{rod}$) acting on each charge. No explicit calculations involving torques or centers of mass, etc. That all emerges as an implicit result of the central force actions. You just explicitly need all those central forces, easy and straightforward to model. Moreover, you can now model not only dipoles, but dipoles plus "free-floating" point charges (i.e., not all your $q$'s need be connected by edges). That would be like an ionized plasma with dipoles in it.
So it's significantly easier to program this way, and lots more general with respect to the kinds of physical situations you can model. Sounds like win-win to me:)
Edit: Extended reply to op's comment below...
Very fundamentally, note that real-life dipoles don't have "virtual rods" constraining the motion of their charges. Consider, for example, a polar water molecule. There are no "rods" -- just forces such that the appropriate molecular Hamiltonian leads to a dipole after solving Schrodinger's equation. Ditto any other kind of dipole -- no "rods". And the center of mass is always an empty point in space, i.e., no physical particle's there. So there's never any force actually acting on the center of mass, simply because there's nothing there to be acted on.
Equations describing center-of-mass motion are themselves simply derived from first principles, i.e., in this case from $\vec F=m\vec a$ and $F=kq_1q_2/r^2$. Generally speaking, if you start from first principles, you can symbolically derive center-of-mass motion and everything else. Therefore, if you numerically model from first principles, your results will necessarily and automatically agree with those symbolically-derived equations.
Consider the off-topic example of Newtonian gravity. Suppose you model the motion of two masses $M\gg m$, say Sun and Earth. If you just use the usual $F=GMm/r^2$, you can forget all about Kepler's laws. The Earth will automatically trace an ellipse around the Sun (with the Sun at a focus), the Earth's areal velocity will automatically be a constant, and ditto $T^2\sim R^3$. You don't have to explicitly build Kepler's laws into your modelling program -- they just emerge from the first-principles dynamics. Indeed, they were derived from the first-principles dynamics (well, okay, Kepler came up with them before Newton, but you get the point:).
So whenever possible, numerical modelling's best approached from first principles. Had you been modelling, say, a gas with $\sim$Avogadro's# particles, then you'd have to model the "emergent" thermodynamic behavior. But I take it your ensemble's small enough for a first-principles approach.
