\documentstyle[12pt]{article}
%\documentstyle[makeidx,psfig,12pt]{article}
\setlength{\topmargin}{-0.90in}
\setlength{\topskip}{0.0in}
\setlength{\topsep}{0.0in}
\setlength{\topsep}{1\topsep}
\setlength{\partopsep}{1\partopsep}
\setlength{\headsep}{0.25in}
\setlength{\headheight}{0.5in}
\setlength{\textheight}{8.90in}
\setlength{\oddsidemargin}{0.5in}
\setlength{\evensidemargin}{0.5in}
\setlength{\textwidth}{5.5 in}
\setlength{\textfloatsep}{1\textfloatsep}

\def\textfraction{.1}
\def\topfraction{.9}

%newcommand for index
%\newcommand{\nn}[1]{#1n}

% My command for a horizontal line
\newcommand{\myrule}{\rule{\textwidth}{0.5mm}}

%\makeindex

\pagestyle{headings}
\begin{document}
\author{by \\
        \ \\ 
        Ziyad S. Hakura \\
        Jim Leathrum \\
        Bill Elliott \\
        Christophe G. Lambert \\
        Bill T. Rankin \\
        and \\
        John A. Board, Jr. \\
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\ 
        \ \\
       Copyright \copyright\ 1995 Duke University \\ All Rights Reserved}
\title {Parallel Multipole Treecode Algorithm \\ PMTA 4.1 \\ Manual}
\date{}
\maketitle
\thispagestyle{empty}
\newpage
\begin{center}
\Large
{Usage Restrictions}
\normalsize
\end{center}

 This code is not in the public domain.  It may, however, be freely
distributed and incorporated into other programs provided:
\begin{enumerate}
\item the copyright notice remains intact, 
\item proper attribution is given to the authors, 
\item no commercial use is made of the program, and 
\item the documentation (PMTAdoc.ps, README\_pmta) accompanies any 
    derivative program. 
\end{enumerate}
\newpage


\section{Introduction}

    This description assumes familiarity with multipole-accelerated
algorithms for the electrostatic N-body problem.
This program computes the Coulomb forces and potentials for a 
system of charged particles; this version is optimized for reasonably
uniform spatial distributions of particles.  The program can also be used
for gravitational problems.  PMTA 4.0 is an entirely new program which
replaces Duke's earlier fast multipole method codes.  PMTA 4.0 is
significantly faster and much more flexible than our earlier codes;
as a new program, this release should be treated as a test release.

    The program is designed to efficiently execute various 
multipole-accelerated algorithms for the N-body problem.  
This particular release can execute
three different algorithms: the Fast Multipole Algorithm (FMA), an
$O(n)$ algorithm 
developed by Greengard and Rokhlin, the Multipole Treecode Algorithm (MTA),
an $O(n \log(n))$ method developed at Duke, and the 
Enhanced Multipole Treecode Algorithm (Enhanced MTA), an $O(n)$ method
which is also developed at Duke.  MTA is a generalization of the
tree code developed by Barnes and Hut (BH).  MTA computes interaction
lists for all cells at the lowest level of spatial decomposition using a
Barnes-Hut-like Multipole Acceptance Criterion (MAC).  In this new MTA 
algorithm, the upward pass is identical to the upward pass of the FMA.  
The downward pass of MTA differs from the FMA in that it only computes forces 
based on the interaction lists of the cells at the lowest level of spatial
decomposition.  It differs from BH in that interactions are on a cell-to-cell
basis and not on a particle-to-cell basis.  
The interaction lists ensure that each particle
interacts, either directly or through a multipole transformation,
with all other particles in the space.  Enhanced MTA is a more 
elaborate version of MTA.  This algorithm allows cells at higher levels
of decomposition to interact with distant cells thus saving the
child cells from interacting individually with the distant cells.  This 
idea is borrowed from the FMA and allows for the linear time compexity of
the two algorithms.

    This code also provides the option of manipulating multipole expansions
in Fourier space; the FFT-enhanced multipole methods (FMA, MTA or Enhanced MTA)
run faster than the original real-space methods at no cost in accuracy but
with significantly higher memory requirements.

Version 4.1 adds periodic boundary conditions, allowing the user to
compute forces and potentials due to macroscopic shapes formed of
replicated unit cubes.  For a given $k$, the potential and forces due
to a $3^k \times 3^k \times 3^k$ region of unit cells can be computed.
Also, computations can be done on macroscopic spheres and ellipsoids, and within
near-infinite square, circular, and elliptical surfaces formed of unit
cubes replicated along two of the three axes.  


    If there are any questions concerning this program, please 
contact Dr. John Board (jab@ee.duke.edu) at the Department of Electrical
Engineering, Duke University.


\section{Compiling}


    The PMTA code can be complied as a stand-alone program for experimenting
with N-body algorithms (a simple main program is provided for this purpose),
or it can be incorporated into a library for invocation from other codes,
i.e. MD programs.

    PMTA can be compiled as a serial program on a workstation or 
as a parallel program on a Kendall Square Research KSR1 shared-
memory computer.

\subsection{Serial Stand-alone compilation}
    To compile the stand-alone version for serial execution, type
\begin{verbatim}
    make 
or  make serial
\end{verbatim}
We recommend the gnu C compiler (gcc) for serial versions.  The program
is known to compile and run correctly with cc on RS-6000's and gcc
on DEC Alpha workstations.  Please report problems with specific platforms
and/or compilers to the authors.  The code is known to raise a number of
ANSI warnings at this point.

\subsection{Parallel Stand-alone compilation}
To compile the stand-alone version for parallel execution on the KSR-1, type
\begin{verbatim}
    make parallel
\end{verbatim}
In either case, the resulting executable is called pmta.  The parallel
version does not support periodic boundary conditions.  There are
plans to rewrite the parallel shared-memory code in a more portable
manner.


\subsection{Library version}
        PMTAmain.c demonstrates the setup and calling conventions of
PMTA, and could be suitably modified to include PMTA as an electrostatics subroutine of a
larger program such as a Molecular Dynamics package. 

        To create a PMTA library, compile all ``.c'' files except
PMTAmain.c, and archive them into a library which will be linked with
your own object files.  Then to incorporate PMTA
routines in your own code you must:
\begin{enumerate}
\item \#include ``PMTAlib.h'' into any files that call PMTA routines.

\item Set up a particle table by calling either 
GenerateParticles(), or ReadParticles(), or with a routine of your own.
GenerateParticles() allocates a particle table, and generates random
particles with random charges within the cubelength specified.
ReadParticles() reads particles from the PDB formatted file using the
given file descriptor.

If you want to just assign positions and charges yourself, do
the following:
\begin{verbatim}
      ParticleTable = (PTablePtr) malloc(NumParticles * sizeof(PTable));      

      for(i=0; i< NumParticles; i++){
         ParticleTable[i].p.x =    
         ParticleTable[i].p.y =    <---- set to appropriate values
         ParticleTable[i].p.z =
         ParticleTable[i].p.q =
      }
\end{verbatim}
\item Call InitializeCoulomb() with the parameters you want.  The
parameters are described in more detail later. 

\item Call Coulomb(ParticleTable)

\item Take your forces answers from the ParticleTable, update the positions
      and goto step 4.

\end{enumerate}


\section{Command line interface}

The command line interface to the supplied main program (PMTAmain.c)
allows the user to evaluate the Coulomb (or Newtonian) interaction
within a fixed or random system of particles.  For a complete 
listing of command line options, type
\begin{verbatim}
    pmta
\end{verbatim}
since without any arguments, the program recognizes that the user 
wishes to see a help screen.

The arguments are listed below with a brief description, and in 
particular cases the range of acceptable values.

Note that this version of PMTA does not read in the values of the
particle charges from the pdb format.

The following arguments are common to both the parallel and serial 
version:

\begin{verbatim}
OPTION    PARAMETERS              DESCRIPTION AND RANGE

One the of following two arguments is necessary but both cannot be
used in the same simulation run:

-a        Theta (double)          Specifies the Multipole Treecode
                                  Algorithm.  Theta is a parameter
                                  used by the Barnes-Hut-like methods
                                  to compute the interaction
                                  lists for all cells at the lowest
                                  level of decomposition (see paper by B+H).
                                  This parameter is specific to the B+H-like
                                  methods and does not appear in the FMA.
                                  Suggested range: 0.5 - 4.0.
                                  (Our definition of theta is slightly different
                                  from the original B+H definition for
                                  efficiency, but the same range of values is
                                  appropriate.)

-pfma     NONE                    Specifies the Fast Multipole 
                                  Algorithm.

Each of the following three arguments is necessary:

-l        Length of Cube          The length of the cube, centered
          (double)                at zero, that contains all the
                                  particles of interest.  If 
                                  particles are read in from a file,
                                  the value entered must be scaled
                                  to the same units as the particle
                                  positions.  Any number
                                  greater than zero is acceptable.
                                  Suggested value: 1.0

-t        Number of Tree Levels   Number of levels in the oct tree.
          (integer)               The first level is level 0.  For
                                  some number x which is input, 
                                  the levels 0,1,2,...,(x-1) will
                                  be used.  Any number greater than
                                  one is acceptable.
                                  Suggested range: depends on number of
                                  particles; practical maximum of about 6.
                                  (The number of cells at a particular
                                  level, w, is 8^w.  Hence if 5 levels
                                  are used, there will be 8^(5-1) = 4096
                                  cells at the lowest level of 
                                  decomposition.)

-p        Number of Multipole     Specifies the number of multipole
          Terms (integer)         expansion terms.  Any number 
                                  greater than zero is acceptable for
                                  real-space algorithms.
                                  For FFT-enhanced methods, only 
                                  4, 8, 12, and 16 may be used.
                                  Suggested range: 4 - 16

The following arguments are optional:

-h        NONE                    Use the Enhanced MTA algorithm.  The
                                  [ -a ] option above must also be
                                  specified together with a Theta
                                  parameter.  A similar Theta parameter
                                  as the MTA aglgorithm is recommended.

-f        NONE                    Use the FFT subroutines in
                                  computing the multipole-to-local 
                                  transformations.  Applies to FMA, MTA and
                                  Enhanced MTA methods.  
                                  This options restricts the 
                                  acceptable values of multipole terms
                                  to 4, 8, 12, and 16.

-n        Number of Particles     Generate a uniform distribution of
          (integer) and           random particles with charges set
          Random seed (integer)   randomly to either +1 or -1.
                                  If this option is not specified then
                                  the program, by default, will read
                                  particle information in pdb format
                                  from standard input.  It will 
                                  randomly assign charges of +1 and -1
                                  to these particles.

-d        NONE                    Compute the exact forces and
                                  potentials for all particles using
                                  the direct method, O(N^2).  The
                                  program performs an error analysis
                                  as a result of this option.  Is also
                                  compatible with -k option, but not
                                  -s or -ellipse.
                                  
-r        NONE                    Print the results for the first
                                  twenty particles.

-o        Percentage Relative     Print detailed information on
          Potential Error         particles which are outliers.
          (double) and            Particles that exceed either
          Percentage Relative     threshold are considered outliers.
          Force Error             The [ -d ] oprion must also be
          (double)                specified on the command line.

-v        NONE                    Print verbose activity messages.

-k        Size of macroscopic     New option, allows periodic 
          region 3^k X 3^k X 3^k  macroscopic region of unit
          (k an  integer)         cells to be computed.  Either
                                  -pfma or -h options must be
                                  selected for this to work. Only
                                  works with serial version.

-s        Radius of macroscopic   Used in conjunction with -k option.
          sphere/disk.            k must be chosen to enclose the
                                  macroscopic sphere.  This option
                                  is slower than straight cubic
                                  macroscopic boundaries, and is
                                  really only here for experimental
                                  purposes.  If code compiled with
                                  -DSURFACE, the unit cube will only
                                  be replicated in the x and y
                                  directions in the shape of a disk 1
                                  unit cube thick.

-ellipse <rx> <ry> <rz>           Dimensions of macroscopic ellipsoid.
                                  Similar to the -s option, only have
                                  a macroscopic ellipse.  Needs -k
                                  with macroscopic cube big enough to
                                  enclose largest ellipse
                                  dimension. If code compiled with
                                  -DSURFACE, the unit cube will only
                                  be replicated in the x and y
                                  directions in the shape of a
                                  elliptical disk 1 unit cube thick.
                                
    
The -s and -ellipse options have not been extensively tested, so be
careful.  For numerical stability when using the FFT [ -f ] option, the
distance units used must be selected such that the length of the
cube containing the simulation is roughly 1E-2 to 1E2 distance units.
In other words, using meters instead of Angstroms as distance units
for a simulation of molecules will not work for the FFT mode of
computation (although it does work for the conventional method of
computation).


The following t ariuments are available only for the parallel version
for the KSR-1 and are optional:

-b        NONE                    Use the Load Balanced parallel model
                                  of execution.  If this option is
                                  not specified then the original
                                  Spatial Decomposition model is 
                                  applied.  The Load Balanced model
                                  can execute on any number of 
                                  processors.  Two features of this
                                  model are static load balancing
                                  and no redundant computations.

-e        Number of processors    Specifies that the program should
          (integer)               not execute on any greater than
                                  this number of processors.  The
                                  program only runs on a power of
                                  two number of processors except
                                  if the Load Balanced parallel model
                                  is executed.
                                  If this option is not specified, 
                                  then, by default, the program 
                                  utilizes the largest power of two 
                                  number of processors within the 
                                  processor set or the total number
                                  of processors within the processor
                                  set if the Load Balanced model is
                                  selected.

-m        NONE                    Print Pmon performance analysis for
                                  a single timestep.

The program prompts the user for any parameters that are required and
missing from the command line.

Example command lines:

pmta -a 1.5 -l 1.0 -t 4 -p 8 -n 30000 1 -d -r >test1

pmta -a 1.5 -l 1.0 -t 4 -p 8 -n 30000 1 -d -r -f >test2

pmta -pfma  -l 1.0 -t 4 -p 8 -n 30000 1 -d -r >test3

pmta -pfma  -l 1.0 -t 5 -p 8 -n 30000 1 -d -r -f -k 2 >test4

\end{verbatim}

\section{Constants}


    The forces and potentials are computed without any constants and
will need to be scaled to the particular units required in a given
application.


\section{Output}

\subsection{Cell Interaction Analysis}

    For cells at the lowest level of decomposition, information is
printed on the number of double direct and single direct interactions
per cell.  Double direct interactions are those direct interactions 
that occur between the particles of two cells that reside on the same 
processor.  Double direct interactions are only recorded on one of the 
two cells' interaction lists.  The double direct interaction list count 
also includes the interaction of a cell with itself.
    Single direct interactions are those direct interactions that occur 
between the particles of two cells that reside on different processors.
Single direct interactions only update the forces and potentials for the
particles that reside locally.


\subsection{Error Analysis}

    Mean relative errors are computed by dividing the mean absolute 
error by the mean absolute value.  Maximum relative errors are computed 
by dividing the worst case absolute error by the mean absolute value.


\subsection{Timing Analysis}

    The first three timing categories listed in the Timing Analysis, 
are not part of a timestep.  All the remaining categories listed 
makeup a single timestep and start with the "Reading Particle Positions"
category.


\section{Outline of the Procedural Interface}

    To call these routines from another program, a procedural interface
is used.  The procedural interface consists of four calls: InitializeCoulomb()
Coulomb(), CoulombDirectBalance() and CoulombMultipoleBalance(). 
InitializeCoulomb() must be called first to specify the parameters of the 
simulation.  For each timestep, Coulomb() is called to compute the 
electrostatic force and potential for all particles.  These calls are 
described in detail in the next two sections.  The CoulombDirectBalance()
and CoulombMultipoleBalance() apply only to the parallel model.  Presently,
CoulombMultipoleBalance() is called once during initialization when
the Load Balanced model is selected and the CoulombDirectBalance() 
procedure is called once at the beginning of every timestep.
Both calls are described in greater detail in the last two sections.
They are not expected to be used externally unless the code is modified.

    For the procedural interface, the three files PMTAmain.c, PMTAfile.c
and PMTAanalysis.c may be omitted from the compilation.
NOTE: the file PMTAparallel.c MUST be included even if the program
is compiled for serial execution.

    Although no analysis is performed on the simulation for this
interface, the timers nevertheless operate in this mode.  
Whenever Coulomb() is called, the per timestep timers are reset to zero
at the beginning of the timestep.  Thus timing data can be accessed 
externally through the Timing structure as follows:

\begin{verbatim}
(double) Timing[int processor-id][int category number][int user or wall]

processor-id    integer between 0 and (the number of processors - 1).

category number integer between 0 and 11 - the numbers correspond 
                to the categories listed in the Timing Analysis output 
                and have the same ordering.
                (NOTE that these numbers correspond to different times
                than in version 3.0.  The PMTA.h file contains a set
                of definitions that have not changed and denote the
                timing data of the various procedures.)

user or wall    0 specifies user time, 1 specifies wall time.


InitializeCoulomb()
===================

NOTE that this call has been modified since versions 3.0 and 4.0.

void InitializeCoulomb(double cubelength, int treelevels, int fourier,
    int pterms, int algorithm, double theta, int processors,
    int direct, int particles, int balance, int verbose, int pmon,
    PTablePtr particletable, int k, int SphereRadius, 
    int EllipseX, int EllipseY, int EllipseZ){

This function should be called with the same values that would be
used in the command-line interface.

ARGUMENT      EQUIVALENT COMMAND-       VALID RANGE OF VALUES
              LINE ARGUMENT

cubelength    -l <Length of Cube>       must be greater than 0, and
                                        enclose all particles when cube
                                        centered at 0,0,0.

treelevels    -t <No. of Tree Levels>   must be greater than 1. 3
                                        gives levels 0,1,2 for 8^2 leaves

fourier       -f                        0 do not use FFT subroutines.
                                        1 use FFT subroutines.

pterms        -p <No. of Terms>         must be greater than 0.
                                        if fourier is 1, must be one
                                        of 4, 8, 12, or 16. 

algorithm     -pfma                     0 Multipole Treecode Algorithm.
                                        1 Fast Multipole Algorithm.
                                        2 Enhanced MTA.

theta         -a <Theta>                B&H Theta parameter. If algorithm
                                        is 0 or 2, must be greater than zero.

processors    -e <No. of Processors>    if less than or equal to zero, 
                                        then use the default processor 
                                        set size.

direct        -d                        0 do not compute exact
                                        particle forces and potentials.
                                        1 compute exact results with O(n^2)
                                        algorithm, in addition to computing the 
                                        approximate values.

particles     -n <No. of Particles>     must be greater than zero.

balance       -b                        0 use Spatial Decomposition 
                                        parallel model.
                                        1 use Load Balanced parallel model.

verbose       -v                        0 be quiet. 1 report on what 
                                        program is doing.

pmon          -m                        Pmon performance: 0=no 1=yes 
                                        (ignored for serial execution)

particletable                           This structure must contain 
                                        the initial positions of the
                                        paricles in the first
                                        timestep. 

k             -k                        Size of macroscopic periodic
                                        simulation (3^k X 3^k X 3^k)

SphereRadius  -s                        -1.0 for no macroscopic sphere,
                                        otherwise, sphere radius.

EllipseX      -ellipse                  -1.0 for no macroscopic ellipse,
EllipseY                                else x,y,z semiaxis dimensions
EllipseZ                                of ellipse.


Example call to InitializeCoulomb():

InitializeCoulomb(1.0,4,1,8,2,1.5,1,0,30000,0,1,0,MyTable,2,-1.0,-1.0,-1.0,-1.0);

This corresponds to a serial simulation of 30,000 particles using the
Multipole Treecode Algorithm and the FFT subroutines, with eight
multipole terms, four tree levels, and a cubelength of one,
periodically replicated in a 9X9X9 volume of unit cells.  MyTable is
of type PTablePtr as described in the next section.

Coulomb()
=========

typedef struct vector {
    double  x,y,z;
} Vector;

typedef struct ptable {
    Vector p;                   /* position */
    int c;                      /* cell identifier */
    double q;                   /* charge */
    Vector fdirect;             /* exact force by direct method */
    double vdirect;             /* exact potential by direct method */
    Vector f;                   /* force */
    double v;                   /* potential */
} PTable;

typedef PTable *PTablePtr;

void Coulomb(PTablePtr particletable);

particletable is a pointer to a linear array of PTable structures 
whose length is equal to the number of particles in the simulation.

Coulomb() uses the position and charge information to update the
fields c, f, and v.  If it was specified in the call to 
InitializeCoulomb() that the results be computed exactly (direct is 1), 
then the fields fdirect and vdirect will also be updated.

CoulombMultipoleBalance()
=========================

void CoulombMultipoleBalance(void);

Called by InitializeCoulomb().  This function implements static
load balancing of the M2L transformations.  This code does not
optimize M2L load balancing based on the presence of empty cells.


CoulombDirectBalance()
======================

void CoulombDirectBalance(void);

Called by Coulomb().  This function implements static
load balancing of the direct particle-particle computations.  
This code considers the positions of the individual particles.
Hence, it can even load balance reasonably non-uniform particle 
distributions.

\end{verbatim}
\section{Reporting Problems}

Please e-mail any questions, problems, bug reports, suggestions, etc.
about this program to: 

\begin{verbatim}
pmta@ee.duke.edu
\end{verbatim}

For more information about multipole algorithm research from the Scientific 
Computing group at the Duke University EE Department, visit our Web page:

\begin{verbatim}
http://www.ee.duke.edu/Research/SciComp/SciComp.html
\end{verbatim}

Available at the Web site are source code files and technical reports explaining 
various aspects of multipole algorithms.
%\bibliographystyle{alpha}
%\bibliography{master}

\newpage

%\printindex
\end{document}

