Primary author: Hans Johnston (johnston@math.umass.edu) Department of Mathematics University of Massachusetts, Amherst This file is the partial documentation for TREECODE3D_TARG, a Fortran90 subroutine for approximating the total electrostatic energy potential and forces of N interacting charged source particles at M target 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_TARG TREECODE3D_TARG 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) and partial charge q_i (i=1,...,N) at the target points Y_k=(x_k,y_k,z_k) (k=1,...,M). For each k the potential at Y_k is given by the scalar N ----- \ q_i PENG(Y_k) = | ----------- (1), / |X_i - Y_k| ----- i=1 where |X_i - Y_k| denotes the Euclidean distance between X_i and Y_k. The electrostatic force on the kth particle is given by the vector F(Y_k)=(f_x,f_y,f_z), and defined by N ----- \ (X_i - Y_k) F(Y_k) = | q_i * ----------- (2), / |X_i - Y_k|^3 ----- i=1 Direct application of either formula (1) or (2) has a computational complexity proportional to O(N^2). TREECODE3D_TARG 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_targ.f : Driver program for testing the treecode subroutine treecode3d_targ.f. treecode3d_targ.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_targ.f treecode3d_targ.f AMD PC running Linux with Portland Group compilers : -------------------------------------------------- pgf90 -fast -byteswapio treedriver_targ.f treecode3d_targ.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_targ.f below). The potential energy or force is approximated by the treecode3d_targ 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_targ.f calling statement: --------------------------------------------------------------------- treedriver_targ.f : -------------- The program treedriver will prompt the user for the following: NUMPARS : (INTEGER) Number of particles THETA : (REAL) The BMAX MAC acceptance parameter. (For more details of the BMAX MAC see the description of THETA below in the section on treecode3d_targ.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 NUMTARS : (INTEGER) For the test program we choose a subset comprised of the first NUMTARS particles to approximate (1) and/or (2). DIST_TOL : (REAL >= 0.0) Distance tolerance used to exclude pairwise interactions unless the distance between them is greater than DIST_TOL. For example, DIST_TOL=0.0 for the test example provided so particle self-interaction is excluded since the targets particles are a subset of the source particles. NOTE: IF THE INTERESTED USER IS SURE THAT THE TARGET AND SOURCE PARTICLE SETS DON'T INTERSECT, THEN THE DISTANCE CHECKS SHOULD BE REMOVED FOR GREATER EFFICIENCY!! treecode3d_targ.f : ------------- The calling statement for the treecode3d_targ 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_TARG(xtar,ytar,ztar,numtars,dist_tol,pengtar, & x,y,z,q,numpars,fx,fy,fz,iflag,order,theta, & shrink,maxparnode,timetree,prinfo,farrdim,numpars) Input: XTAR,YTAR,ZTAR : (REAL*8 one-dimensional arrays of length NUMTARS) Cartesian coordinates of the target particles NUMTARS : (INTEGER) number of target particles DIST_TOL : (REAL*8) Potential energy/force contribution is excluded if any particle lies within DIST_TOL of a target particle. (NOTE: this parameter is only used in the case of a direct computation) X,Y,Z : (REAL*8 one-dimensional arrays of length ARRDIM) Cartesian coordinates of the aource particles. Q : (REAL*8 one-dimensional array of length ARRDIM) Corresponding charge of each particle. NUMPARS : (INTEGER) number of source 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. (IF IFLAG=1 then FARRDIM is ignored, otherwise it should be set to NUMTARS) ARRDIM : (INTEGER) declared length of arrays X,Y,Z, and Q Output: PENGTAR : (REAL*8 one-dimensional arrays of length FARRDIM) k-th entry is potential energy at k-th target particle. FX,FY,FZ : (REAL*8 one-dimensional arrays of length FARRDIM) force on particle is given by (FX(k),FY(k),FZ(k)) for k=1,...,NUMTARS. TREETIME : (REAL*8) time (in seconds) for the tree algorithm computation.