Primary author: Hans Johnston (johnston@math.umass.edu) Department of Mathematics University of Massachusetts, Amherst This file is the partial documentation for TREECODE3D, a Fortran90 subroutine for approximating the total electrostatic energy potential and forces of N mutually interacting charged particles using an adaptive treecode. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% README file for TREECODE3D TREECODE3D is a Fortran90 subroutine for approximating the electrostatic energy potential and force of N mutually interacting charged particles with Cartesian coordinates X_i=(x_i,y_i,z_i) (i=1,N) and partial charges q_i (i=1,N). The potential energy PENG is given by the scalar N-1 N ----- ----- \ \ q_i*q_j PENG = | | ----------- (1), / / |X_i - X_j| ----- ----- i=1 j=i+1 where |X_i - X_j| denotes the Euclidean distance between X_i and X_j. The electrostatic force on the ith particle is given by the vector F(X_i)=(f_x,f_y,f_z), and defined by N ----- \ (X_i - X_j) F(X_i) = q_i | q_j * ----------- (2) , / |X_i - X_j|^3 ----- j=1 j ~= i Direct application of either formula (1) or (2) has a computational complexity proportional to O(N^2). TREECODE3D evaluates the potential energy (1) and electrostatic force (2) using a hierarchical oct-tree algorithm with body-cell forces approximated via Taylor expansions. The Taylor coefficients are computed by a recurrence relation. This approach reduces the computational complexity to O(NlogN). NOTE: Please include the following references in any work that utilizes this code: (1) Duan, Z.-H., Krasny, R.: An adaptive treecode for computing nonbonded potential energy in classical molecular systems. J. Comput. Chem. {\bf 22} (2001) 184-195 (2) Lindsay, K., Krasny, R.: A particle method and adaptive treecode for vortex sheet motion in 3-D flow. J. Comput. Phys. {\bf 172} (2001) 879-907 (3) Lindsay, K.: A three-dimensional Cartesian tree-code and applications to vortex sheet roll-up. Ph.D. Thesis, University of Michigan (1997) Summary of files in ./example and ./aux : ----------------------------------------- ./example treedriver.f : Driver program for testing the treecode subroutine treecode3d.f. treecode3d.f : Treecode subroutine. xyx3d.dat : Unformatted big endian binary file containing the coordinates of 50,000 particles uniformly distributed in [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5] q3d.dat : Unformatted big endian binary file containing 50,000 partial charges uniformly distributed in [-0.5,0.5]. input_50kp : Input data file provided for testing the treecode for computing the potential energy. The files xyz3d.dat and q3d.dat provide the particle data. output_50kp : Program output using input_50kp. These results were produced using the the SUNWspro f90 compiler on a 750MHZ SparcIII Sun Blade 1000 with 512MB of memory. input_50kpf : Input data file provided for testing the treecode for computing the potential energy AND force. The files xyz3d.dat and q3d.dat provide the particle data. output_50kpf : Program output using input_50kpf. These results were produced using the the SUNWspro f90 compiler on a 750MHZ SparcIII Sun Blade 1000 with 512MB of memory. ./aux gendata3d.m : MATLAB m-file script provided to generate randomly distributed coordinate and charge data for an arbitrary number of particles. The user is prompted for N, the number of particles, and the code generates both an xyz3d.dat and q3d.dat file. The supplied files xyz3d.dat and q3d.dat (in ./example) are sufficient to test the code for any N <= 50,000. Compiling and Testing the code ------------------------------ Compiling : The source code is written in Fortran90. Two compile line examples are provided below, which both produce an executable a.out. SUN BLADE : --------- f90 -fast -O4 treedriver.f treecode3d.f AMD PC running Linux with Portland Group compilers : -------------------------------------------------- pgf90 -fast -byteswapio treedriver.f treecode3d.f The second example has been provided to emphasize that the 64-bit REAL binary data files provided in ./example are big endian format. Thus on some architectures a compiler flag (e.g. -byteswapio for the pgf90 compiler) may be needed to ensure proper conversion to your machine's native 64-bit REAL type. Testing : After compiling, to test the code type either of the following: (potential energy ONLY) ./a.out < input_50kp (potential energy AND force) ./a.out < input_50kpf This executes the code with N=50,000 and 5 other user specified parameters (see the description of inputs for treedriver.f below). The potential energy or potential energy force is approximated by the treecode3d subroutine and then computed exactly by direct summation (formulas (1) and (2)). Both the absolute and relative error in the infinity norm of the potential energy and the three force components (viewed as N vectors) are output, as well as the total time (in seconds) for each of the calculations. To ensure that the code is functioning properly on your machine compare your screen output with the output in the file output_50kp or output_50kp. Description of treedriver.f inputs and treecode3d.f calling statement: --------------------------------------------------------------------- treedriver.f : -------------- The program treedriver will prompt the user for the following: NUMPARS : (INTEGER) Number of particles THETA : (REAL*8) The BMAX MAC acceptance parameter. (For more details of the BMAX MAC see the description of THETA below in the section on treecode3d.f.) ORDER : (INTEGER) Order the Taylor expansions used for the approximation if the MAC is accepted. MAXPARNODE : (INTEGER) When the tree is constructed if the number of particles in the current node exceeds MAXPARNODE then the box defined by the node particles is subdivided. A smaller MAXPARNODE generally results in a deeper tree structure. SHRINK : (INTEGER) Adaptive switch used in the oct-tree construction. If equal to 1 then upon subdivision of a node the node bounds are taken the smallest Cartesian box containing the particles in the node. Otherwise each box is divided into 8 sub-boxes of equal area as determined by the midpoints in each coordinate direction. IFLAG : (INTEGER) Flag determining what is to be computed: IFLAG=1 : potential energy IFLAG=2 : potential energy AND forces treecode3d.f : ------------- The calling statement for the TREECODE3D subroutine is as follows: !!!!!!!! IMPORTANT: If you only need the potential energy set IFLAG=1. You MUST still pass the FX, FY, and FZ variables. However, they can be scalars in which case you should set FARRDIM=1. CALL TREECODE3D(x,y,z,q,numpars,fx,fy,fz,peng,iflag,order,theta, & shrink,maxparnode,timetree,prinfo,farrdim,arrdim) Input: X,Y,Z : (REAL*8 one-dimensional arrays of length ARRDIM) Cartesian coordinates of the particles. Q : (REAL one-dimensional array of length ARRDIM) Corresponding charge of each particle. NUMPARS : (INTEGER) number of particles IFLAG : (INTEGER) flag indicating what is to be computed: 1 - potential energy 2 - potential energy and force ORDER : (INTEGER) order of Taylor expansion used in the body-cell approximation. THETA : (REAL*8) opening angle used in the BMAX MAC (multi-pole acceptance criteria). If P is the position of the target particle and R is the "radius" (distance from the geometric center C of the box to its furthest corner) of the current box under consideration then the Taylor expansion is used if R/|P-C| < THETA. SHRINK : (INTEGER) adaptive switch used in the oct-tree construction. If equal to 1 then upon subdivision of a box the box bounds are taken the smallest Cartesian box containing the particles in the node. Otherwise each box is divided into 8 sub-boxes of equal area. MAXPARNODE : (INTEGER) When the tree is constructed if the number of particles in the current node exceeds MAXPARNODE then the box defined by the node particles is subdivided. A smaller MAXPARNODE generally results in a deeper tree structure. PRINFO : (INTEGER) switch for outputting information in this routine. If set to 1 then information is written to stdout. FARRDIM : (INTEGER) declared length of arrays FX, FY, and FZ. ARRDIM : (INTEGER) declared length of arrays X,Y,Z, and Q. Output: PENG : (REAL*8) potential energy FX,FY,FZ : (REAL*8 one-dimensional arrays of length FARRDIM) force on ith particle is given by (FX(i),FY(i),FZ(i)). TREETIME : (REAL*8) time (in seconds) for the tree algorithm computation.