\documentclass{nic-series}

\begin{document} 

% all words in headings should be capitalised except for conjunctions, prepositions, definite and indefinite articles 
% unless they appear at the beginning
\title{Molecular Dynamics Algorithms for \\Massively Parallel Computers}

% list of authors and institutes
\author{Norbert Attig \inst{1} \and 
        Marius Lewerenz \inst{2} \and \\
	Godehard Sutmann \inst{1} \and 
        Reiner Vogelsang \inst{3}}

% list of authors for the table of contents, please use abbreviations for the first names
\authortoc{N. Attig, M. Lewerenz, G. Sutmann, R. Vogelsang}

% list of organisations, numbered automatically in the order of your input
\institute{Central Institute for Applied Mathematics,\\
         Forschungszentrum J\"ulich, 52425 J\"ulich, Germany\\
         \email{\{n.attig, g.sutmann\}@fz-juelich.de}
         \and
         LADIR/Spectrochimie Mol\'eculaire, Universit\'e Pierre et Marie Curie, \\
         4, Place Jussieu, 75252 Paris C\'edex 05, France \\
         \email{lewerenz@spmol.jussieu.fr}
         \and
         Silicon Graphics GmbH, Professional Services,\\
         Am Hochacker 3, 85630 Grasbrunn, Germany \\
         \email{reiner@munich.sgi.com}
          }

\maketitle

\begin{abstracts}
A new modular program system, DMMD (Distributed Memory Molecular
Dynamics) has been designed which employs modern
elements of the Fortran~90 language for abstraction and encapsulation.
The program uses the message passing paradigm which is applicable
to both distributed and shared memory systems.
Communication is handled in a special interface library which
efficiently translates high level Fortran~90 concepts into the more
primitive addressing mechanisms of current message passing protocols.
Furthermore it allows transparent switching between message passing 
libraries. This makes it a powerful, flexible tool for parallel Fortran 90
applications. The modular concept of DMMD allows easy modification and
extensions of its capabilities towards new interaction models or
integrators. The design has been tested in a serial version which also
serves as a reference point for the parallel version which will be
used to benchmark several strategies for data and work distribution.
We report first results obtained with DMMD on a Cray T3E system
for a simple Lennard-Jones mixture and comparisons with other
massively parallel molecular dynamics codes.
\end{abstracts}

\section{Introduction}

Within the John von Neumann Institute for Computing (NIC)
the Central Institute for Applied Mathematics (ZAM) acts as a national
supercomputer centre and provides access to high performance computer
resources and programming expertise to a wide scientific community.
\newline

$\vdots$
\newline

These applications will be outlined in Sec.~\ref{attig_sec_motive}.
Strategies for the parallel evaluation of forces are discussed in
Sec.~\ref{attig_sec_forces}.
The design concept of the DMMD molecular dynamics program is described
in Sec.~\ref{attig_sec_design} and the key features of a general Fortran~90
communication interface and other software engineering considerations
are the subject of Sec.~\ref{attig_sec_software}.
A description of important features of the main hardware platform, 
CRAY T3E, and benchmark results are given in Secs.~\ref{attig_sec_hardware} 
and \ref{attig_sec_benchmark}, respectively.

% to ensure the clearness of references of all contributions the name of the author should be attached to the labels
\section{Motivation}
\label{attig_sec_motive}
One challenge in modern high performance computing is the development
of programs which are able to simulate systems on a microscopic
basis, i.e. atomic resolution, on {\em large} time and spatial scales, 
where {\em large} means micrometers and microseconds. 
This will open the field to treat problems like crystal growth, crack
propagation or protein folding which could, until now, only be described
with various simplifications. 
Fig.~\ref{attig_systems_figure} gives a schematic overview over different 
simulation methods and their range of applicability.
\newline
        
$\vdots$

\begin{figure}[t]
\begin{center}
%  [width=...] will enlarge or reduce the figure based on the value for width
%  [width=\textwidth] will use the full textwidth for the figure
\includegraphics[width=8cm]{systems.pdf}
\caption{\label{attig_systems_figure}
Schematic picture of the application domains of 
quantum molecular dynamics (QMD), classical molecular dynamics (CMD), 
Brownian dynamics (BD) and hydrodynamics on time and spatial scales.
}
\end{center}
\end{figure}

\section{Parallel Force Evaluation in Molecular Dynamics}
\label{attig_sec_forces}
This section describes different types of algorithms for short
range interactions in parallel molecular dynamics calculations.
The two types of algorithms which we discuss in detail are the
particle decomposition (PD) and spatial decomposition (SD) methods.
The characteristics of these algorithms are in principle very
different and they have their own advantages depending on the
problem which is to be solved.

$\vdots$

\subsection{Neighbour Lists}
Each MD step requires the calculation of the forces acting on a
particle due to all other particles in the simulation box.
\newline

$\vdots$
\subsection{Particle Decomposition Method}
{\bf PD1:}
In this algorithm each processor stores the data of the coordinates,
velocities and forces of a subset of particles $N_p = N/P$,
where $N$ is the total number of particles of the system and
$P$ is the number of processors.
It is assumed here and in the following that $\mbox{mod}(N/P)=0$ but
it is easy to generalise the algorithm to non-homogeneous distributions.
In order to test for mutual interactions between particles the
coordinates of all particles, stored on other processors,
are communicated to the local PE.
The most natural way to realise the communication would be to
perform on each processor P-1 send/receive operations
so that the coordinates of all particles in the system are
duplicated on each PE.
A more efficient way which needs only $\log_2(P)$ send/receive 
operations, was proposed by Fox et al.\cite{attig_fox88} and is illustrated
in Fig.~\ref{attig_pd1_figure}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{pd1.pdf}
\caption{\label{attig_pd1_figure}
Example for the PD1 algorithm on 8 processors. 
Left: The force loop requires the whole vector ${\bf x}$ on all 
processors. 
The equations of motion are solved on one processor with a 
subset ${\bf x}_i$, ${\bf v}_i$, ${\bf F}_i$ on processor
PE$_i$. 
Right: The procedure to copy the whole vector ${\bf x}$ on one PE
requires only $\log_2(P)$ steps.
}
\end{center}
\end{figure}
$\vdots$
\newline
\noindent
With this convention one may subdivide the force loop into
loops of equal length through
\newline                                                         

\noindent
{\tt
                offset = 0                                   \newline
                do ij = 1,n\_loop                            \newline
\hspace*{0.2cm}   offset = offset + 1                        \newline
\hspace*{0.2cm}   do i = 1,N$_p$                             \newline
\hspace*{0.4cm}     j = i + offset                           \newline
\hspace*{0.4cm}     $\delta {\bf x} = {\bf x}_i - {\bf x}_j$ \newline
\hspace*{0.6cm}     $\vdots \qquad$ ({\it force evaluation}) \newline
\hspace*{0.2cm}   enddo                                      \newline
                enddo                                        \newline
}

\noindent
where $n_{loop} = (N_p -1)/2$. The performance of PD1 is summarised 
in Tab.~\ref{attig_Tab1}.

\begin{table}[ht]
\centerline{
\begin{tabular}{ll}
\hline \\ [-1.2ex]
1. All-to-all communication of coordinates: \hspace*{0.8cm} &   $N-N_p$  \\
   \hspace*{0.35cm}
   local copy of total vector ${\bf x}$ on each  PE    &                 \\ [0.4ex]
2. Construction of neighbour lists:                    &                 \\
   \hspace*{0.35cm}
   {\tt if( $L<4R_c$)} Verlet list                     &   $N_p(N -1)$   \\
   \hspace*{0.35cm}
   {\tt if( $L\ge 4R_c$)} linked cell list             &   $N_p$         \\ [0.4ex]
3. Calculation of force vector for $i \in N_p$         &   $N_p$         \\ [0.4ex]
4. Update particle positions                           &   $N_p$         \\ [0.4ex]
\hline
\end{tabular}
}
\caption{\label{attig_Tab1}
  Scaling behaviour with number of particles for PD1.}
\end{table}


\section{Molecular Dynamics Program Design}
\label{attig_sec_design}

\section{Software Engineering}
\label{attig_sec_software}

\section{Hardware Considerations}
\label{attig_sec_hardware}
The DMMD program runs on workstations (SUN, IBM R50, DEC) as well
as on powerful massively parallel processor (MPP) systems like 
CRAY~T3E and Origin~2000.
Program development is done on workstations, while full production 
runs with large numbers of particles are performed on MPP systems only;
benchmarks and test simulations run on all platforms available.
Target platforms for DMMD are primarily distributed memory MPP systems,
e.g. the CRAY~T3E, which have the following advantages over
shared memory systems:
(i) programs have the most general structure, i.e.
algorithms and methods developed for distributed MPP systems are also
applicable to shared memory machines, which is not true {\em vice versa}!
(ii) the CRAY~T3E is currently one of the most powerful supercomputers
available\footnote{Two large T3E systems are installed at the
Forschungszentrum in J\"ulich.
The first one is equipped with 512 processors of type T3E-600 and a
main memory of 128 MByte per PE.
The second one consists of 512 processors of type  T3E-1200E
with 512 MByte main memory per PE.},
and it is especially suited to tackle large-scale MD problems.
This was demonstrated by {\em world record simulations} of the
groups of Trebin\cite{attig_Trebin} and of Kollman\cite{attig_Kollman}.
\newline

$\vdots$

\section{Benchmarks}
\label{attig_sec_benchmark} 

\subsection{Proper MD Programs for Comparison}
\label{attig_sim_md_progs}
The simplest benchmark programs available for our purpose are the ones
from Steve Plimpton which are freely available on the 
web\cite{attig_plimpton_web}.
They simulate Lennard-Jones mixtures with different parallelisation
strategies: particle distribution (PD1, PD2), domain (SD) and force
decomposition\cite{attig_plimpton95}.
\newpage
LAMMPS (developed at Sandia National Laboratory within the CRADA 
project\cite{attig_LAMMPS_web}), a pure MD program, partially contains
these algorithms. It is designed to model
large molecular systems and it is also described on the web.
\newline

$\vdots$

\subsection{The Benchmark Run}
\label{attig_bench_run}
Our benchmark problem is a simulation of a simple Lennard-Jones fluid
with the parameters given in Tab.~\ref{attig_Tab2}. 
\newline

$\vdots$

\begin{table}[h]
\centerline{
\begin{tabular}{p{0.2cm}ll}
\hline \\ [-1.2ex]
$\bullet$  & \multicolumn{2}{l} {System of short range interacting particles:} \\
           \hspace*{0.35cm}
           & \multicolumn{2}{l} {Argon-Krypton mixed Lennard-Jones (LJ) fluid} \\ [0.4ex]
$\bullet$  & \multicolumn{2}{l} {NVE ensemble}                                 \\ [0.4ex]
$\bullet$  & \multicolumn{2}{l} {Global physical parameters:}                  \\ [0.4ex]
           & density:         & $n = 0.0181\: \text{\AA}^{-3}$                 \\ [0.4ex]
           & temperature:     & $T = 116\: \text{K}$                           \\ [0.4ex]
           & LJ parameters:   & $\sigma_{\text{Ar}}   = 3.405\: \text{\AA}$,
                                $\epsilon_{\text{Ar}} = 119.8\: \text{K}$      \\
           \hspace*{0.35cm}
           &                  & $\sigma_{\text{Kr}}   = 3.633\: \text{\AA}$,
                                $\epsilon_{\text{Kr}} = 167.0\: \text{K}$      \\ [0.4ex]
           & mass:            & $m_{\text{Ar}} = 39.948\: \text{amu}$,
                                $m_{\text{Kr}} = 83.8\: \text{amu}$            \\ [0.4ex]
           & cutoff radius:   & $r_{\text{cutoff}} = 2.5 \: \sigma$            \\ [0.4ex]
           & skin radius:     & $r_{\text{skin}}   = 0.5 \: \sigma$            \\ [0.4ex]
           & time step:       & $\Delta t = 20\: \text{fs}$                    \\ [0.4ex]
$\bullet$  & \multicolumn{2}{l} {Global simulation parameters:}                \\ [0.4ex]
           & system:          & CRAY T3E-1200, streams enabled                 \\ [0.4ex]
\hline
\end{tabular}
}
\caption{\label{attig_Tab2} Benchmark Parameters}
\end{table}

\subsection{Further Benchmark Remarks}
\label{attig_bench_remarks}
In a real production run neighbour lists are obligatory.
The list update has only to be performed in certain time intervals, 
depending on the skin radius and the diffusivity of the system.
DMMD does not yet contain neighbour lists but we have tested this 
technique with Plimptons' kernels where additionally Newton's 
third law is not applied (algorithm PD1) and neighbour lists are 
updated at regular intervals of 20 time steps.
Tab.~\ref{attig_Tab3} reports the results of these runs.

\begin{table}[ht]
\centerline{
\begin{tabular}{|r|r|r|r|r|r|r|r|r|}
\hline
     & \multicolumn{4}{c|} {N = 2048} & \multicolumn{4}{c|} {N = 97556} \\\hline
     & \multicolumn{2}{c|} {particle} &
       \multicolumn{2}{c|} {domain}   &
       \multicolumn{2}{c|} {particle} &
       \multicolumn{2}{c|} {domain}   \\[-1.0ex]
     & \multicolumn{2}{c|} {decomposition} &
       \multicolumn{2}{c|} {decomposition} &
       \multicolumn{2}{c|} {decomposition} &
       \multicolumn{2}{c|} {decomposition} \\\hline
     & Verlet & linked & Verlet & linked & Verlet & linked & Verlet & linked \\[-1.0ex]
     &        & \multicolumn{1}{c|} {cell} & & \multicolumn{1}{c|} {cell} &
              & \multicolumn{1}{c|} {cell} & & \multicolumn{1}{c|} {cell} \\\hline
\#PEs& \multicolumn{4}{c|} {time/step/N [$\mu$s]} & \multicolumn{4}{c|} {time/step/N [$\mu$s]} \\\hline
   1 & 49.31 & 41.71 & 20.65 & 17.63 & 927.9 & 39.10 & 761.2 & 16.94 \\\hline
   2 & 24.81 & 21.01 & 12.74 & 10.30 & 463.3 & 19.54 & 245.4 &  8.87 \\\hline
   4 & 12.67 & 10.76 &  7.34 &  6.01 & 231.9 &  9.88 &  78.8 &  4.68 \\\hline
   8 &  6.57 &  5.61 &  4.19 &  3.56 & 116.1 &  5.09 &  25.5 &  2.44 \\\hline
  16 &  3.55 &  3.07 &  2.21 &  2.10 &  56.7 &  2.66 &   8.0 &  1.28 \\\hline
  32 &  2.09 &  1.86 &  1.25 &  1.29 &  29.3 &  1.51 &   2.6 &  0.66 \\\hline
  64 &  1.34 &  1.20 &  0.82 &  0.82 &  14.8 &  0.86 &   0.8 &  0.34 \\\hline
 128 &  1.00 &  0.94 &  0.62 &  0.64 &   7.5 &  0.55 &   0.3 &  0.19 \\\hline
\end{tabular}
}
\caption {\label{attig_Tab3} Plimpton-Benchmarks with 2048 and 97556 particles}
\end{table}

First of all, the table clearly shows that domain decomposition is 
faster than particle distribution, nearly independent of the number 
of particles.
\newline

$\vdots$
\section{Concluding Remarks}
We have described a molecular dynamics package for systems with short
range interactions. 
The program is written in a strong typing Fortran 90
dialect and designed in a strictly modular form.
A message passing wrapper has been developed which supports different 
message passing libraries, e.g. MPI, and avoids common problems 
between Fortran 90 and communication libraries. 
Several types of algorithms for short range
interactions were discussed and the performance of the program was
compared to Plimpton's benchmark kernels\cite{attig_plimpton95}.
\newline

$\vdots$

\section*{Acknowledgements}
We are indebted to G. Groten for many fruitful discussions about the
intricacies of Fortran~90 and for the coding of a significant part
of the interface library.
The benchmark computations were performed with a grant of computer time
provided by the VSR of the Forschungszentrum J\"ulich.

%  attach the name of the author to the bibitems to ensure the clearness of references of all contribution
\begin{thebibliography}{10}

\bibitem{attig_plimpton95}
S.~J.~Plimpton,
{\em Fast parallel algorithms for short-range molecular dynamics},
\Journal{J.~Comp. Phys.}{117}{1--19}{1995}.

\bibitem{attig_fox88}
G.~C.~Fox, M.~A.~Johnson, G.~A.~Lyzenga, S.~W.~Otto, J.~K.~Salmon, and D.~W.~Walker,
{\em Solving Problems on Concurrent Processors: {V}olume 1}
(Prentice Hall, Englewood Cliffs, NJ, 1988).

\bibitem{attig_Trebin}
MD simulation of an FCC lattice with 1.213.857.792 atoms, see
\url{http://www.itap.physik.uni-stuttgart.de/$\sim$joerg/imd.html}

\bibitem{attig_Kollman}
MD simulation of a 1 $\mu$s trajectory of a small protein in water:
(a) Y.~Duan, L.~Wang and P.~A.~Kollman,
\Journal{Proc. Natl. Acad. Sci. USA}{9}{9897}{1998}.
(b) Y.~Duan and P.~A.~Kollman,
\Journal{Science}{282}{740}{1998}.

\bibitem{attig_plimpton_web}
\url{http://www.cs.sandia.gov/$\thicksim$sjplimp/main.html}

\bibitem{attig_LAMMPS_web}
\url{http://www.cs.sandia.gov/$\thicksim$sjplimp/crada.html}

\end{thebibliography}

\end{document} 
