How best to simulate n-body systems in a functional way?

Question

I want to simulate a N-body system in the most functional way possible. The system I'm modeling are metallic colloidal particles in a solution of peroxide water, as here http://arxiv.org/pdf/1502.07144.pdf. Due to the small size of the particles, the system is in a low Reynold's regime, so there is no inertia. Furthermore, the surface of the particles react with the fluid and there is creation and consumption of hydrogen (H) on their surface. Gradients on the concentration of H will produce forces on the particles, which depend on their material parameters. The particles can show self-organization and motion.

In general, the long range force between the particles is given by the combined action of the hydrodynamic forces and the gradient of the H concentration field, which is yet to be experimentally characterised. We assume that the force between the particles is a radially symmetric, pair-wise function that decays as $1/r^2$. In principle this force can be modified by changing the chemistry of the system and/or the composition of the particles, say, by adding salts and changing the potential to a Yukawa potential instead of $1/r^2$. I'm interested on finding the structures formed by these colloidal particles that show dynamic function, e.g., oscillators, by modifying the forces or the material parameters.

The systems I study have generally only a few particles, so I'm not so concerned with an efficient implementation of the force computation. Rather I want a very abstract code that allows me to eventually run several thousands of simulations of small systems with different parameters, analyse the structures found, and re-run with slightly different parameters —- think of genetic algorithms. A function that runs fast would be nice though.

This is one example of the structures found for two kind of particles interacting:

I already have a C code that integrates the system and works well, but I find the problem of solving a general N-body problem in Mathematica using its built-in integrators interesting in own right and as an example of how straight-forward such code might be in a functions style. (I took some ideas and code from 38687). Further, the Mathematica code is 100 times shorter (albeit slower) and could be easily expanded to simulate molecules, granular matter or the like by simply changing the intereaction force and leaving all the technical details to a fast implementation to compute the force between the particles. I think this would be great for teaching molecular dynamic simmulations for example.

I've been using Mathematica for some years for data analysis but recently started to study more in depth the functional way of programming, so excuse me if the code is not elegant, it's working though. This is where I am at the moment.

Define some particles, define a certain functional form for the force and a interaction matrix that is not necessarily symmetric. (If you simulate gravity, for example, you have a symmetric matrix, in diffusophoretic systems you generally don't.)

Then integrate with NDSolve.

Clear[force, r, type]
numbodies = 9
pos0 = RandomReal[{-10, 10}, {numbodies, 2}]

type = 1 + Boole[OddQ[#]] & /@ 
  Range[numbodies] (*type of particle, we start with 2 but they could \
be more, keepp it balanced to form structures*)

matrixB = 
 2.51 {{1, -1}, {-.55, 
    1.1}} (*interaction matrix between the particles*)

force[j_] := 
 Sum[(Normalize[
     r[j][t] - r[i][t]]) (matrixB[[type[[i]], type[[j]]]]/
      EuclideanDistance[r[j][t], r[i][t]]^2 +
     1/EuclideanDistance[r[j][t], r[i][t]]^12), {i, 
   Delete[Range[numbodies], 
    j]}];(*diffusophoretic force and leonard jones repulsion*)

odesys = Table[{r[j]'[t] == force[j], r[j][0] == pos0[[j]]}, {j, 
   numbodies}];
depvars = Flatten[Table[{r[j]}, {j, numbodies}]];
tfin = 500;
sol = NDSolve[odesys, depvars, {t, 0, tfin}][[1]];

pos = Array[r, {numbodies}] /. sol;
Animate[Show[
  ParametricPlot[#[t] & /@ pos, {t, Max[0, tfin1 - 50], tfin1}, 
   ImageSize -> 400, Frame -> True, PlotRange -> {-20, 20}], 
  Graphics[
   MapIndexed[{Hue[.35 type[[First[#2]]] ], Disk[#1[tfin1], 0.5]} &, 
    pos]]], {tfin1, 0.1, tfin, 1}]

So, my questions are:

Is this a proper way to functionally code an N-body integration? Are there better/more-elegant ways of doing it? Can the performance be improved? Is the Animate of the plot re-computing the solution for each final time it shows? Can I pass parameters to the interaction matrix without Mathematica complaining? My objective is to be able to answer questions concerning the interaction matrix and the structures formed, so this code is just a little part of big multi-dimensional minimization, say, to find the matrix/force law/initial-condition that maximizes the rotational energy in the system.

(if you are wondering why the values of the interaction matrix, they allow for the formation of a pure oscillator in some cases.)

Feel free to edit the title of the question if it's not well chosen. I learnt a lot while writing this question so maybe it's not the best phrasing now that I have NDSolve working.

EDIT

I think I haven't been clear enough with the generality of the problem. This is basically the same code to solve a 3-body gravitational problem.

Clear[force, r, type, v]
numbodies = 3
pos0 = RandomReal[{-10, 10}, {numbodies, 2}]
vel0 = RandomReal[{-.1, .1}, {numbodies, 2}]
vel0 = (# - Mean[vel0]) & /@ vel0
type = 1;
matrix = {{1, 1}, {1, 1}} (*interaction matrix between the particles*)
force[j_] := 
 Sum[-(Normalize[r[j][t] - r[i][t]]) matrix[[type, type]]/(.1 + 
      EuclideanDistance[r[j][t], r[i][t]]^2), {i, 
   Delete[Range[numbodies], 
    j]}];(*gravitational force with a kernel*)

odesys = Table[{v[j]'[t] == force[j],
   r[j]'[t] == v[j][t],
   r[j][0] == pos0[[j]],
   v[j][0] == vel0[[j]]}, {j, numbodies}];


depvars = Flatten[Table[{r[j], v[j]}, {j, numbodies}]];

tfin = 200;
sol = NDSolve[odesys, depvars, {t, 0, tfin}][[1]];

pos = Array[r, {numbodies}] /. sol;


Animate[Show[
  ParametricPlot[#[t] & /@ pos, {t, Max[0, tfin1 - 50], tfin1}, 
   ImageSize -> 400, Frame -> True, PlotRange -> {-20, 20}], 
  Graphics[MapIndexed[{Hue[.35], Disk[#1[tfin1], 0.5]} &, 
    pos]]], {tfin1, 0.1, tfin, 1}]

And this is how it looks (don't know why the flinchs appear when I export to gif...):

In this case the interaction matrix is just a scalar, and it's where the Gm1m2 should go. Also the particles have only one type. The equation of motion in this case is a bit more complicated since you need equations for both position and velocity, but the way of solving the problem is exactly the same. I hope this clarify a bit my question.

Answer 1

In version 12.0 you can use NBodySimulation for that:

data  = NBodySimulation[
   "InverseSquare", {<|"Mass" -> 1, "Position" -> pos0[[1]], 
     "Velocity" -> vel0[[1]]|>,
    <|"Mass" -> 1, "Position" -> pos0[[2]], 
     "Velocity" -> vel0[[2]]|>,
    <|"Mass" -> 1, "Position" -> pos0[[3]], 
     "Velocity" -> vel0[[3]]|>}, tfin];

Visualize the data:

ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, tfin}]