\chapter{THE MOLECULAR DYNAMICS PROBLEM}

\label{chapter_md}

  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.75\textwidth]{chap_md/md.png}
    \caption{Particles with Momentum and Weak Attraction}
  \end{centering}
  \end{figure}

\section{Introduction}

A great deal is known about the structure and behavior of atoms and
molecules, as long as we consider them one at a time.  Some properties of
a material can be explained directly by corresponding features that can be
noticed in a single atom.  But there are other phenomena, such as phase
transitions, the formation of surfaces, the effect of catalysts, and various
mixing effects, whose understanding can only arise from considering
large groups of molecules having a typical range of energies, and interacting
by collisions and forces of attraction or repulsion.  

The field of molecular dynamics provides a simulated laboratory for such studies.
Instead of atoms or molecules, we consider a collection of abstract {\it{particles}},
which don't react chemically.  We think of them as point masses which move
though space, driven by their momentum.  A pair of particles exert a mutual
repulsive force when close, and a weaker attractive force when further apart.

A computation might involve hundreds of thousands of particles.  The strategy
requires creating a series of ``snapshots'' of the particle positions over time.
However, to achieve accurate results over a time interval of interest, a huge
number of snapshots must be taken.  To create the next snapshot requires 
selecting each particle, measuring its distance to all other particles,
summing up the resultant forces, and then using Newton's law to move the
particle a tiny distance.  The same process must be carried out for each of
the particles before the snapshot can be completed.  If we have $n$ particles
to study, then we have roughly $n^2$ pairs of molecular interactions to sum
up on a snapshot.  Thus, a molecular dynamics code can provide a problem with
quadratic growth; as $n$ increases, our problem quickly grows from difficult
to intractable.  (For this very reason, sophisticated molecular dynamics codes 
cut down the work by concentrating on the local molecular interactions, while
approximating the weaker effect of far away molecules.)

Packages for doing heavy duty molecular dynamics simulation include AMBER,
CHARMM, DL\_POLY, Gromacs, and NAMD.

\section{The Calculation}

Suppose that we have $n$ particles, and that at some starting time $t_0$ 
we are given, for each particle $i$, its mass $m_i$, and the initial 
position, velocity, and acceleration vectors $p_i(t), v_i(t), a_i(t)$.  
Our goal is to model the behavior of the particles
up to some final time $t*$ by measuring the forces on the particles, and allowing
them to respond over a short time step.  If our steps are small enough and
our approximation scheme accurate, then we can get reliable estimates of the
behavior of the particle system over the time range of interest.

In order
to determine what is going to happen at time step $t+dt$, we need to compute
the force vector $f$.  In the molecular dynamics model, the
force on particle $i$ is the sum of the individual forces exerted by all
the other particles in the system.  The force $f_{i,j}$ exerted on particle $i$ by
particle $j$ depends only on the distance $d_{i,j}=||p_i-p_j||$ between the two particles,
and its value is $f_{i,j} = V'(d_{i,j})$, that is, the derivative of the potential
function.  This step is carried out in the function {\it{compute\_new()}}.

The potential function used in the {\it{md()}} program is the sine-squared
potential $V(d)=\sin^2(\min(d,\frac{\pi}{2}))$, which is completely zero beyond 
a distance of $\frac{\pi}{2}$, and attractive at closer range.

  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.49\textwidth]{chap_md/sinsq_potential.png}
    \includegraphics[width=0.49\textwidth]{chap_md/sinsq_derivative.png}
    \caption{The $\sin^2$ Potential V(d) and Derivative V'(d)}
  \end{centering}
  \end{figure}

(A more realistic potential function is known as the Lennard-Jones potential.  
Particles repel each other when close because of the Pauli principle, and 
attract at a distance because of the van der Waals force.)

Once the force vector $f$ has been computed, it is possible to determine
the position, velocity and acceleration vectors for snapshot at time $t+dt$
using the values known at time $t$.  One technique for doing this is known 
as the {\it{Velocity Verlet}} algorithm:
\begin{align*}
  p(t+dt) =& p(t) + v(t) dt + \frac{1}{2} dt^2 \\
  a(t+dt) =& f(t) / m \\
  v(t+dt) =& v(t) + \frac{1}{2} ( a(t) + a(t+dt) )
\end{align*}
This step is carried out in the function {\it{update()}}.

Once the new position, velocity and acceleration have been computed,
a new snapshot is complete.  The program can gather statistics from this
new data, such as the pressure, temperature, kinetic and potential energy,
before beginning the computation of the next step.

\section{Estimating the Work}

A typical run of our particular program might involve 1,000 particles for 
500 time steps.  This seems to mean that a million computations of some sort
will be involved, which suggests our program is not going to run in the blink
of an eye!  In fact, if we are a little more careful, we can make a reasonable
estimate of how the work and time will grow with the values of $np$, the number
of particles, and $step\_num$ the number of steps.  

Assume we are using $np$ particles, and taking $step_num$ time steps.
The basic calculation in the code computes the total force on one particle $p_i$
at the $k$-th time.  This requires summing up the individual force from each particle
$p_j$.  If it costs $c$ units of work to calculate an individual force, then it
will cost $(np-1)*c$ units to get the total force on particle $p_i$, and thus
$np*(np-1)*c \approx np^2*c$ units to compute the total force on all the particles, 
at the $k$-th time.  To do this for each time step will cost us a total work
of about $step\_num*np^2$ units.  This means that the problem is linear in
time steps, but {\it{quadratic}} in the number of particles.  You should be
warned that problem work depends quadratically on a parameter, that the work 
can very quickly explode beyond your computer's capacity.

If we really believe that the work is linear in time, then we can estimate the
work involved in a 1,000 particle and 500 timestep run by using {\bf{tic}}
and {\bf{toc}} to time 1 step of a 1,000 particle run, and multiplying by 500.  
While we're at it, we can do this experiment for a range of particle sizes as well.

\begin{tabular}{|r|r|r|r|}
\hline
NP & 1 step time & 500 step estimate & 1,000 step estimate \\
\hline
1,000 &  1s &  8m & 17m\\
2,000 &  4s & 33m & 1h 6m\\
4,000 & 14s & 1h 57m & 3h 53m\\
8,000 & 58s & 8h 3m & 16h 6m\\
16,000 & 3m 51s & 1d 8h & 2d 16h\\
\hline
\end{tabular}

Our estimates are a warning that this program can eat up a lot of time,
especially if we try to increase the particle size.  This suggests that it would
be worth our trouble to try to rewrite the code in a way that speeds it up.
It's hard to be certain about what to try; before changing our code, the next
thing we should do is find out what lines of code are being used the most.  This
is where our work is being done, and this is where we should concentrate our efforts.

{\it{Typical calculations with a production molecular dynamics code might involve 
20,000 to 3,000,000 particles and take 10,000 to 50,000 time steps!  The problem of
controlling the computational cost becomes so serious that special algorithms have
been developed to cut down on the quadratic cost of the particle count parameter.}} 

\section{Profiling the Sequential Code}

The molecular dynamics calculation that we are using as an example is unlike
most of the other example programs we have encountered so far.  It's longer,
made up of multiple functions, and involves loops that look complicated.
That's why it is not easy to point to any particular line of the program and
assume that speeding up that line will speed up the entire program.  We often
aren't able to form a realistic assessment of the importance of the individual
portions of our code.

In practice, however, it is often the case that, when a large problem is being
solve, most of the computational work occurs in just a few program locations; 
If we can locate those locations, and try to improve their operation, we can achieve
a huge improvement in performance.  And, conversely, it is often the case that
programmers waste time optimizing a portion of their program that contributes almost
nothing to the overall cost.   Since we only
have a certain amount of time to spend modifying and testing the code, it
is wise to identify those sections containing most of the work, and then
concentrate on them alone.

Thus, we propose to run the program sequentially, using a realistic problem
size, and try to spot the busiest loops, where optimization might help.  Fortunately,
MATLAB comes with a very useful tool for this purpose, the {\it{program profiler}}.
To use the program profiler for the {\bf{md}} program, we simply type the following
three commands:
\begin{lstlisting}
profile on
md
profile viewer
\end{lstlisting}
When the program is complete, the command {\bf{profile viewer}} opens up a
report on the program which lists the functions that were called (both user
and internal MATLAB functions), the number of times called, the total time
spent in the function or its subfunctions, and the ``self time'' spent precisely
in that function itself.  

Now it's normal for the program to take somewhat longer to run if the profiler is
on, because 
When the profiler is turned on, MATLAB is doing more work that just executing your program;
it also keeps track of which statement is begin carried out, and how long it took.  Typically,
this adds a noticeable amount of time to the original running time.  However, there's
something about the {\bf{md}} program that makes the profiled version much, much
slower.  Running with 1,000 particles and 100 steps, the program took about 90 seconds to
run, but with the profiler turned on, the program used more than 5,000 seconds.  This is
surprising, but it's actually an indirect indication that there's something wrong with
the program.  We won't try to explain the time discrepancy further, and will concentrate
on the report itself.
%
%  Labels for figures must go inside the caption field or LATEX doesn't
%  set them correctly.
%
  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.99\textwidth]{chap_md/md_profile1.png}
    \caption{A portion of the profile report for {\bf{md}} \label{fig-md-profile1} }
  \end{centering}
  \end{figure}
  
Figure~\ref{fig-md-profile1} displays the initial part of the profile report for the {\bf{md}} program,
and indeed {\bf{md}} shows up on the first line, called once (when we invoked it
from the command line), and with a total time equal to the the program running time.
However, the self time is negligible - that is, as soon as we call {\bf{md}}, 
{\bf{md}} is calling other functions, and that is where the work is really going on.

Moving to the second line of the report, we see the item {\bf{md$>$compute}}, meaning
that we are reporting on time spent in the {\bf{compute}} function, called by {\bf{md}}.
This is called 101 times (once for start up, and once for each of 100 steps).
The total time is again very large, but now notice that the self time is about the
same.  In other words, we've found where the work is actually being carried out.
Although {\bf{md}} includes functions {\bf{initialize()}} and {\bf{update()}}, 
the report indicates that their contribution to the total cost is negligible.

Clearly, our concern should focus on {bf{compute()}}.  Fortunately, the profiler
can give us more detailed information of what is going on in that function.  We simply
click on the {\bf{md$>$compute}} line of the report, and we are transferred to
a document that analyzes the {\bf{compute()}} function in greater detail (Figure~\ref{fig-md-profile2}).

  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.99\textwidth]{chap_md/md_profile2.png}
    \caption{A portion of the profile report for {\bf{md$>$compute}} \label{fig-md-profile2} }
  \end{centering}
  \end{figure}
  
Now we see individual lines of the code labeled with the number of times they
were called, and their contributions to the overall time.  From this report,
it is clear that 50\% of the total program time is spent executing 5 lines of code
(really just 3 lines, since the other two lines are {\bf{end}} statements).
\begin{lstlisting}
  f = zeros ( nd, np );
  pot = 0.0;
  for i = 1 : np
    for j = 1 : np
      if ( i ~= j )
        for k = 1 : nd
          rij(k) = pos(k,i) - pos(k,j);
        end
        d = 0.0;
        for k = 1 : nd
          d = d + rij(k)^2;
        end
        d = sqrt ( d );
        d2 = min ( d, pi / 2.0 );
        pot = pot + 0.5 * sin ( d2 ) * sin ( d2 );
        for k = 1 : nd
          f(k,i) = f(k,i) - rij(k) * sin ( 2.0 * d2 ) / d;
        end
      end
    end
  end
  kin = 0.0;
  for k = 1 : nd
    for j = 1 : np
      kin = kin + vel(k,j)^2;
    end
  end
  kin = 0.5 * mass * kin;
\end{lstlisting}

At this point, we need to get advice from someone familiar with MATLAB optimization
about whether we can improve the function.  However, you should be aware that MATLAB
is not efficient when asked to carry out {\bf{for}} loops that don't actually contain
a lot of work, and are called many times.  If possible, such loops should be replaced
by equivalent vector operations, which MATLAB can execute much more quickly.

On the advice of an experienced MATLAB code, the {\bf{compute()}} function was
revised.  The new function was called {\bf{compute\_new()}}, and it was put inside
a new program called {\bf{md\_new()}}.  The core of the new function is dramatically
changed:

\begin{lstlisting}
  f = zeros ( nd, np );
  pot = 0.0;
  pi2 = pi / 2.0;
  for i = 1 : np
    Ri = pos - repmat ( pos( :, i ), 1, np );
    D = sqrt ( sum ( Ri.^2 ) );
    Ri = Ri( :, ( D > 0.0 ) );
    D = D( D > 0.0 );
    D2 = D .* ( D <= pi2 ) + pi2 * ( D > pi2 );
    pot = pot + 0.5 * sum ( sin ( D2 ).^2 );
    f( :, i) = Ri * ( sin ( 2*D2 ) ./ D )';
  end
  kin = 0.5 * mass * sum ( diag ( vel' * vel ) );
\end{lstlisting}

Now let us redo the simple timing estimate with the revised code to see if the
changes had any effect:

\begin{tabular}{|r|r|r|r|}
\hline
NP & 1 step time & 500 step estimate & 1,000 step estimate \\
\hline
 1,000 & 1/5s &  1m40 & 3m 20s\\
 2,000 &  2/3s & 5m33 & 11m 6s\\
 4,000 & 2s & 16m 40s & 33m 20s \\
 8,000 & 8s &  1h 6m&  2h 13m\\
16,000 & 38s & 5h 16m&  10h 33m\\
\hline
\end{tabular}

After seriously reworking the code, we've managed to get it to run about 4 times faster
on the larger problems.  In one sense, this is wonderful.  In another sense, it's unfortunate
because after working so hard on the {\bf{compute}} function, it's really hard to see how
we can squeeze out any more speed, short of rethinking the algorithm.

\section{Performance of the Parallel Code}

If we want to run bigger problems, or get our results faster, then the natural thing to
consider is to carry out the calculations in parallel.  If we run the profiler on the revised
{\bf{md\_new}}, we will still observe that the {\bf{compute\_new}} function accounts for
almost all the work.  Luckily, the work is contained inside a {\bf{for}} loop that carries
out the force calculation for each particle.  These force calculations are independent.
This means we can set up a parallel version of the program quite simply.
\begin{lstlisting}
  f = zeros ( nd, np );
  pot = 0.0;
  pi2 = pi / 2.0;
  parfor i = 1 : np
    Ri = pos - repmat ( pos( :, i ), 1, np );
    D = sqrt ( sum ( Ri.^2 ) );
    Ri = Ri( :, ( D > 0.0 ) );
    D = D( D > 0.0 );
    D2 = D .* ( D <= pi2 ) + pi2 * ( D > pi2 );
    pot = pot + 0.5 * sum ( sin ( D2 ).^2 );
    f( :, i) = Ri * ( sin ( 2*D2 ) ./ D )';
  end
  kin = 0.5 * mass * sum ( diag ( vel' * vel ) );
\end{lstlisting}

If this code will execute in parallel, then if we decide to run a bigger, longer job,
we 
