# Simulating molecular dynamics efficiently

The following website offers a very nice molecular dynamics simulation: http://physics.weber.edu/schroeder/md/. It is pretty neat and quite a few physical phenomena can be described from that small scale, see the author's flyer. The author mentions that the code can run smoothly in a browser with 1000 particules and a 2000 steps per second (only displaying a few dozens per second, of course). I wanted to simulate something similar with Mathematica but am facing efficiency issues. I am not interested in features such as live manipulation of particules, I mostly want to choose initial conditions, compute the time evolution, and visualize the results later.

This is my code:

n = 200; (* number of particules *)
h = 0.01; (* time step *)
size = 50.; (* size of the box *)

Clear[x, a]
(* initial conditions *)
x[0] = Table[{2 Mod[i, size] - size, Floor[i/50] - size}, {i, n}];
SeedRandom[1234];
v[0] = 0.1*size*RandomReal[{-1, 1}, {n, 2}];
x[1] = x[0] + h*v[0];
(* Verlet integration scheme *)
mem : x[i_] := mem = 2 x[i - 1] - x[i - 2] + h^2*a[i - 1];
mem : a[i_] := mem = Table[forceParticules[x[i], j, epsilon], {j, n}]
+  Table[forceWalls[x[i], j, 100], {j, n}];


Explanations: the particules are initial stacked at the bottom of a square box, with random initial velocities. As in the applet, the time integration scheme is a Verlet integration scheme. The acceleration at time i+1, denoted by a[i], comprises two things:

• the force from surrounding particles (forceParticules). The trick is to consider only the influence of surrounding particles, at distance of at most epsilon, to reduce computation time. The surrounding particles exert of force deriving from a Lennard-Jones potential.
• the force exerted by the boundary of the box (forceWalls).

This is the definition I used, and both are called hundreds of hundreds of times.

r0 = 1; (* diameter of one particle *)
epsilon = 10; (* radius of influence *)
potential[r_] = 4 * ((r0/r)^12 - (r0/r)^6);
dpotential[r_] = potential'[r];
forceParticules[xi_, i_, epsilon_] := Module[{},
close = Nearest[Drop[xi, {i}], xi[[i]], {100, epsilon}];
vecs = # - xi[[i]] & /@ close;
output = Total[dpotential[Norm[#]]*#/Norm[#] & /@ vecs]
]

stiffness[x_] =
Piecewise[{{-100*(x + size), x < -size}, {-100*(x - size),
x > size}}, 0];
forceWalls[xi_, i_, raideur_] := {stiffness[xi[[i, 1]]],
stiffness[xi[[i, 2]]]}


It then suffices to evaluate the x[i] and view the result:

 Table[x[i], {i, 500}];
Animate[ListPlot[x[i], PlotRange -> size*{{-1, 1}, {-1, 1}},
AspectRatio -> 1, PlotStyle -> PointSize[1/size]], {i, 1, 100, 1}]


My problem is that the code runs very slowly (13 seconds for 500 steps with 50 particles). I am miles away from the efficiency of the applet.

• For a start, memoization is the best way to slow it down... Better use arrays instead. Dec 6, 2017 at 18:35
• @HenrikSchumacher Could you explain your claim? From dzone.com/articles/java-8-automatic-memoization, I read that "If it takes a long time to compute the return value, memoization will speed up the program. For functions that take less time to evaluate than fetching the previously stored value from memory, memoization is clearly not a good option._" In my case, computing the value is certainly longer than fetching the stored value. I would think that forceParticle is more responsible than memoization, which is only used 1000 times here. Dec 6, 2017 at 19:01
• Then you take to much time to compute the values... ;o) Memoization in a hash table still has an overhead. To do this is only justified if 1.) it is not predictable at coding time when the value has to be retrieved and 2.) if the values cannot be stored in a predictable ordering. Both of these conditions are not fulfilled: 1.) In each time step, only the values from the preceding time iteration are required. 2.) the time steps have a natural order. Dec 6, 2017 at 19:06
• There are also other issues. Nearest is applied in an inefficient way; in each iteration, a single call to Nearest is required (at most) and the computation of the forces is better delegated to a CompiledFunction. Dec 6, 2017 at 19:08
• @HenrikSchumacher Thank you for the clarification. I'm working on the code accordingly. Regarding your second comment, I don't understand how you can, with a single use of Nearest, find for each particle the position of neighbour particles. Dec 6, 2017 at 19:12

Okay, here is a way to compute the forces much faster: We create a CompiledFunction (called getForces). It eats a list of points in the plane and spits out the net force onto the first point of the list; here the second to last points are supposed to be those points that are so close to the first one that they exert a force onto it.

size = 50.;(*size of the box*)
box = {{-size, size}, {-size, size}};
r0 = 1.;(*diameter of one particle*)

Quiet[Block[{r, x1, x2, y1, y2, xx, yy, force, potential},
xx = {x1, x2};
yy = {y1, y2};
potential = r \[Function] 1/4 ((r0/r)^2 - (r0/r));
force = -D[potential[Sqrt[Dot[xx - yy, xx - yy]]], {xx, 1}];
With[{
f1code = N@force[[1]], f2code = N@force[[2]], slope = 100.,
a1 = N@box[[1, 1]], b1 = N@box[[1, 2]], a2 = N@box[[2, 1]],
b2 = N@box[[2, 2]]
},
getForces = Compile[{{X, _Real, 2}},
Block[{x1, x2, y1, y2, f1, f2},
x1 = CompileGetElement[X, 1, 1];
x2 = CompileGetElement[X, 1, 2];
f1 = slope (Ramp[a1 - x1] - Ramp[x1 - b1]);
f2 = slope (Ramp[a2 - x2] - Ramp[x2 - b2]);
Do[
y1 = CompileGetElement[X, i, 1];
y2 = CompileGetElement[X, i, 2];
f1 += f1code;
f2 += f2code;
, {i, 2, Length[X]}
];
{f1, f2}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
]]];


The following is a (not so pure) function that computes the $n \times 2$ matrix xnew consisting of the new point positions; it uses the matrix x of positions at a given time instance and the matrix xold which represents the positions at the preceeding time instance. xold is handled as side effect which makes it not too nicely, but this way, we can use it in NestList later.

step = x \[Function] (
xnew =  2. x - xold +  h^2 (1./m) getForces[Nearest[x, x, {\[Infinity], epsilon}]];
xold = x;
xnew
);


Setting up the remaining parameters and the initial conditions...

SeedRandom[1234];
n = 200;(*number of particules*)
h = 0.01;(*time step*)
(*initial conditions*)
x0 =   N@Table[{2 Mod[i, size] - size, Floor[i/50] - size}, {i, n}];
v0 = 0.1 size RandomReal[{-1, 1}, {n, 2}];
x1 = x0 + h v0;(*Verlet integration scheme*)
m = ConstantArray[1., n];(*particle masses*)


... and now, we let it run:

timesteps = 10000;
xold = x0;
x = x1;
data = Join[{x0}, NestList[step, x1, timesteps]]; // AbsoluteTiming


{4.58476, Null}

Aha, roughly 2000 iterations per second. That's at least not so much worse than the JavaScript implementation...

And here a visualization:

frames = Map[X \[Function] Graphics[{PointSize[1/size ],
Point[X]}, PlotRange -> box, PlotRangePadding -> 0.1 size
], data];
Export["a.gif", frames[[1 ;; -1 ;; 20]]]


Note that I use a different potential (less singular) in order to get it running. Of course, you can play with the parameters...

Further improvements

In the meantime, I did some hand tuning of the compiled code. The new positions get now computed completely within the CompiledFunction. That's the first new thing. The second is that for each point also all current positions, the preceeding position of this point, and an index list ilist with the indices for the relevant points are handed over. The first index in ilist markes the point for which we want to compute the new position; the other entries mark the points in the region of interaction. This way, we can recycle the index lists ilists that Nearest will produce every now and then in the main loop (see below).

size = 50.;(*size of the box*)
box = {{-size, size}, {-size, size}};
r0 = 1.;(*diameter of one particle*)
Quiet[
Block[{x1, x2, y1, y2, xx, yy, force, potential, r, r2},
xx = {x1, x2};
yy = {y1, y2};
potential = r \[Function] 4 ((r0/r)^12 - (r0/r)^6);
force = Simplify[
-D[potential[Sqrt[Dot[xx - yy, xx - yy]]], {xx, 1}]
/. (x1 - y1)^2 + (x2 - y2)^2 -> r2
] /. Sqrt[r2] -> r;

With[{
f1code = N@force[[1]], f2code = N@force[[2]], slope = 100.,
a1 = N@box[[1, 1]], b1 = N@box[[1, 2]], a2 = N@box[[2, 1]],
b2 = N@box[[2, 2]]
},
getStep =
Compile[{{X, _Real, 2}, {Xold, _Real, 1}, {ilist, _Integer,
1}, {factor, _Real}},
Block[{x1, x2, y1, y2, f1, f2, r, r2, j, i},
j = CompileGetElement[ilist, 1];
x1 = CompileGetElement[X, j, 1];
x2 = CompileGetElement[X, j, 2];
f1 = slope (Ramp[a1 - x1] - Ramp[x1 - b1]);
f2 = slope (Ramp[a2 - x2] - Ramp[x2 - b2]);
Do[
i = CompileGetElement[ilist, k];
y1 = CompileGetElement[X, i, 1];
y2 = CompileGetElement[X, i, 2];
r2 = (x1 - y1)^2 + (x2 - y2)^2;
r = Sqrt[r2];
f1 += f1code;
f2 += f2code;
, {k, 2, Length[ilist]}
];
{
2. x1 - CompileGetElement[Xold, 1] + factor f1,
2. x2 - CompileGetElement[Xold, 2] + factor f2
}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]

]];


Preparing some constants...

SeedRandom[1234];
n = 1000;(*number of particules*)
h = 0.005;(*time step*)

(*initial conditions*)

x0 = DeveloperToPackedArray[
N@Table[{2 Mod[i, size] - size, Floor[i/50] - size}, {i, n}]];
v0 = 0.1 size RandomReal[{-1, 1}, {n, 2}];
x1 = x0 + h v0;

m = ConstantArray[1., n];(*particle masses*)
factors = (h^2./m);
timesteps = 10000;
skip = 10; (*Nearest gets called only every 10th time iteration*)


Main loop

We write directly into a preallocated array (which makes no difference since NestList is cleverly implemented). The major change is that we request only the vertex indices from Nearest (with x -> Automatic). In contrast to the coordinates, these won't change significantly for several iterations. So we have to call Nearest less than once per time iteration. DeveloperToPackedArray seems to be needed because the (listable) getStep is applied to ragged lists (which cannot be packed) and so the output is not packed.

data = ConstantArray[0., {timesteps + 1, n, 2}];
data[[1]] = x0;
data[[2]] = x1;
xold = x0;
x = x1;
ilists = Nearest[x -> Automatic, x, {\[Infinity], epsilon}];
Do[
If[Mod[iter, skip] == 0,
ilists = Nearest[x -> Automatic, x, {\[Infinity], epsilon}];
];
data[[iter]] = xnew = DeveloperToPackedArray[getStep[x, xold, ilists, factors]];
xold = x;
x = xnew;
,
{iter, 3, timesteps + 1}]; // AbsoluteTiming


{4.37231, Null}

The timing is much better than 2000 steps per second now.

I have also prepared a somewhat nicer visualization:

velocities = Sqrt[(Join[{v0}, Differences[data]/h]^2).ConstantArray[1., {2}]];
colorcoords = Rescale[Clip[velocities, {-100, 100}]];
frames = Table[
Graphics[{
PointSize[0.5 r0/size],
Transpose[{ColorData["TemperatureMap"] /@ colorcoords[[i]],
Point /@ data[[i]]}],
}, PlotRange -> box, PlotRangePadding -> 0.1 size,
Background -> Black
],
{i, 1, Length[data], 100}];

Manipulate[frames[[i]],{i, 1, Length[frames], 1}]


• That's really great. I think it should be PointSize[0.5*r0/size] instead of PointSize[r0/size] in the last code block. I'll try to edit your answer with other funny animations if I get some. Thank you very much. Dec 7, 2017 at 18:51
• Also Mod[iter, skip] == 0 instead of == 1, so it works also with skip=1? Dec 7, 2017 at 18:53
• You're welcome! Why not posting examples in an own answer. Answering your own questions is legitimate, you know... Dec 7, 2017 at 18:55
• Yes but nice animations (obtained with your code) are not an answer on how to simulate efficiently. Dec 7, 2017 at 18:57

To expand on @HenrikSchumacher's comment, compare:

r1 = Table[
Nearest[Drop[x[1],{i}], x[1][[i]], {100,10}],
{i,200}
];//AbsoluteTiming

nf = Nearest[x[1]]; //AbsoluteTiming
r2 = nf[x[1], {100, 10}][[All, 2;;]]; //AbsoluteTiming

r1 === r2


{0.004298, Null}

{0.000029, Null}

{0.000463, Null}

True

Also, note that you can have the NearestFunction return distances as well:

r1 = Table[
Norm[#-x[1][[i]]]& /@ Nearest[Drop[x[1],{i}], x[1][[i]], {100,10}],
{i,200}
]; //AbsoluteTiming

nf = Nearest[x[1]->"Distance"];//AbsoluteTiming
r2 = nf[x[1], {100, 10}][[All, 2;;]]; //AbsoluteTiming

r1 === r2


{0.021323, Null}

{0.000049, Null}

{0.000326, Null}

True

• Thank you for the addendum too! Very handy. I don't know how a user can know that given I don't see it the doc, though! Dec 6, 2017 at 20:42
• It's in the documentation for Nearest, e.g., example Dec 6, 2017 at 21:12
• The functionality must have been introduced in 11.1 because I don't see it in 11.0. Dec 6, 2017 at 22:17
• @anderstood Yes, it was. Before 11.1 you have -> Automatic, which is the same as -> "Index"`, but no other property. Dec 7, 2017 at 7:41
• Ah, that explains, why "Index" does not work for me. Thanks for the clarification! Dec 7, 2017 at 18:31