## The Numerical NBody Problem An Introduction to Molecular Dynamics

One way that systems of ODEs arise in the physical sciences is in the description of the motion of N interacting classical particles. Newton solved this problem for the case of two particles (the two-body problem) interacting via a central force. However, for three or more particles there is no general analytic solution, and numerical techniques are of great importance in understanding the motion.

In the numerical method known as molecular dynamics, the coupled equations of motion for the N particles are simply integrated forward in time using Newton's second law for each particle. There is nothing subtle about this—the numerical techniques learned in the previous sections are simply applied on a larger scale. The subtleties only arise when details such as error accumulation, code efficiency, and the like must be considered.

Below, we show how to use Mathematical intrinsic function NDSolve to numerically solve the following N-body problem: For particles at positions r 1(t), r 2(t),..., rN (t), the equations of motion are d2 r. N

where F, is the force between particles i and j, and mi is the mass of particle i.

To complete the problem, we must also specify initial conditions on position and velocity:

The following module MDsol solves this problem numerically:

Cell 1.49

Clear["Global'*"]

Cell 1.50

MDsol[z0_, time_] := Module[{}, npart = Length[z0]/6;

r[i_, t_] = {x[i][t], y[i][t], z[i][t]}; v[i_, t_] = {vx[i][t], vy[i][t], vz[i][t]};

Z[t_] = Flatten[Table[{r[i, t], v[i, t]}, {i, 1, npart}]]; f[t_] = Flatten[Table[{v[i, t],

(Sum[F[i, j, r[i, t]-r[j, t]], {j, 1, i-1}] + Sum[F[i, j, r[i, t]-r[j, t]],

{j, i + 1, npart}]}/mass[[i]]}, {i, 1, npart}]]; ODEs = Flatten[Table[Z"[t][[n]] == f[t][[n]],

{n, 1, 6 * npart}]]; ics = Table[Z[0][[n]] == z0[[n]], {n, 1, 6*npart}]; eqns = Join[ODEs, ics] ;

NDSolve[eqns, Z[t], {t, 0, time }, MaxSteps™ 10*5]]

To understand what this module does, look at the last line. Here we see that NDSolve is used to integrate a list of equations called eqns, that the equations involve a vector of unknown functions Z[t] , and that the equations are integrated from t = 0 to t = time. The definition of Z[t] can be found a few lines higher in the module: it is a list of variables, {r[i,t],v[i,t]}. The ith particle position vector r[i,t] is defined in the second line as having components {x[i][t],y[i][t],z[i][t]}, and the velocity vector v[i,t] has components {vx[i][t],vy[i][t],vz[i][t]}. These functions use a notation we haven't seen before: the notation x[i][t] means the same thing as x[i,t]. The reason we use the former and not the latter notation is due to a vagary of NDSolve: NDSolve likes to work on functions of one variable; otherwise it gets confused and thinks it is solving a PDE. The notation x[i][t] fools NDSolve into thinking of x as a function of a single argument, the time t.

The Flatten function is used in the definition of Z[t] because NDSolve works only on a simple list of unknown functions, without sublists.

The list eqns can be seen to be a concatenation of a list of ODEs called ODEs and a list of initial conditions called ics. The initial conditions are given as an argument to the module, in terms of a list z Q of positions and velocities for each particle of the form z 0 = Flatten[{r ^ V1Q, r 20 , V20,...)]

The Flatten command is included explicitly here to ensure that z0 is a simple list, not an array.

The acceleration function f(t) is defined as and the result is flattened in order for each element to correspond to the proper element in Z. The fact that the sum over individual forces must neglect the self-force term j = i requires us to write the sum as two pieces, one from j = 1 to i - 1, and the other from j = i + 1 to N.

The value of N, called npart in the module (because N is a reserved function name), is determined in terms of the length of the initial condition vector z0 in the first line of the code.

Finally, the module itself is given no internal variables (the internal-variable list is the null set { }), so that we can examine each variable if we wish.

In order to use this code, we must first define a force function E,-(r). Let's

consider the gravitational N-body problem, where the force obeys Newton's law of gravitation:

where G = 6.67 X 10 11 m3/kg s2. We can define this force using the command

F[i_, j_, r_] := -mass[[i]] mass[[j]] r/(r.r) *(3/2)

Here mass is a length-N list of the masses of all particles, and we have set the gravitational force constant G = 1 for simplicity.

Let's apply this molecular dynamics code to the simple problem of two gravitating bodies orbiting around one another. For initial conditions we will choose r 1 = Yj = 0, and r2 = (1,0,0), v2 = (0,0.5,0). Thus, the list of initial conditions is z0 = Flatten[{{0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0.5, 0}}]

Also, we must not forget to assign masses to the two particles. Let's take one mass 3 times the other:

Cell 1.51

Cell 1.52

Cell 1.53

mass

 {{x[1][t] ™ InterpolatingFunction[{ {0., 4-}} , <>] [t]. y[i][t] ™ InterpolatingFunction[{ {0., 4-}} , <>] [t]. z[1][t] ™ InterpolatingFunction[{ {0., 4-}} , <>] [t]. vx[1][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t] vy[1][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t] vz[1][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t] x[2][t] ™ InterpolatingFunction[{ {0-, 4-}} , <>] [t]. y[2][t] ™ InterpolatingFunction[{ {0., 4-}} , <>] [t]. z[2][t] ™ InterpolatingFunction[{ {0., 4-}} , <>] [t]. vx[2][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t] vy[2][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t] vz[2][t] ™InterpolatingFunction[ {{0. , 4.} }. <> ][t]

The result is a list of interpolating functions for each component of position and velocity of the two particles. We can do whatever we wish with these—perform more analysis, make plots, etc. One thing that is fun to do (but is difficult to show in a textbook. is to make an animation of the motion. An example is displayed in the electronic version of the book, plotting the (xy) positions of the two masses at a series of separate times. (The animation can be viewed by selecting the plot cells and choosing Animate Selected Graphics from the Format menu.) Only the command that creates the animation is given in the hard copy, in Cell 1.55.

### Cell 1.55

Table[ListPlot[Table[{x[n][t], y[n][t]}/ .% [[1]], {n, 1, npart}], PlotStyle™ PointSize[0.015], AspectRatio™ 1, PlotRange™{{- .1, 1.2}, {-.1, 1.2}}], {t, 0, 4, .1}];

Another thing one can do (that can be shown in a textbook!) is plot the orbits of the particles in the x-y plane. The parametric plots in Cell 1.56 do this, using the usual trick of turning off intermediate plot displays in order to save space. The mass-1 particle can be seen to move considerably farther in the x-direction than the mass-3 particle, as expected from conservation of momentum. Both particles drift in the y-direction, because the mass-1 particle had an initial y-velocity, which imparts momentum to the center of mass.

Cell 1.56

p1 = ParametericPlot[{x[1][t], y[1][t]}/.S[[1]], {t, 0, 4}, DisplayFunction™Identity];

p2 = ParametericPlot[{x[2][t], y[2][t]}/.S[[1]], {t, 0, 4}, DisplayFunction™Identity,

Show[p1, p2, DisplayFunction™\$DisplayFunction, PlotRange-> {{0, 1}, {0, 1}}, AspectRatio™1];