Primary author: Hans Johnston (johnston@math.umass.edu) Department of Mathematics University of Massachusetts, Amherst Secondary author: Peijun Li (lipeijun@math.purdue.edu) Department of Mathematics Purdue University Copyright (c) 2009. All Rights Reserved. This file is the documentation for TREECODE3D_YUKAWA (v1.1), a Fortran90 subroutine for approximating the total electrostatic energy potential and force of N mutually interacting charged particles undergoing screened Coulomb (Yukawa) interactions 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_YUKAWA (version 1.1 9/15/09) TREECODE3D_YUKAWA is a Fortran90 subroutine for approximating the electrostatic energy potential and force of N mutually interacting charged particles undergoing screened Coulomb (Yukawa) interactions 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 at the ith particle is given by the scalar N ----- \ EXP(-kappa|X_i - X_j|) PENG_i = q_i | q_j ---------------------- , (1) / |X_i - X_j| ----- j=1 j~=i where |X_i - X_j| denotes the Euclidean distance between X_i and X_j, and kappa >= 0 is the Debye-Huckel parameter. The force vector at the ith particle is given by the vector ( FORCE = -gradient(PENG) ) N ----- \ EXP(-kappa|X_i - X_j|)(kappa|X_i-X_j|+1) FORCE_i = q_i | q_j (X_i-X_j) ----------------------------------------, (2) / |X_i - X_j|**3 ----- j=1 j~=i Direct application of (1) or (2) has a computational complexity proportional to O(N^2). TREECODE3D_YUKAWA evaluates the potential energy (1) and force (2) using an adaptive hierarchical oct-tree algorithm with body-cell interactions approximated via Taylor expansions. The Taylor coefficients are computed by a recurrence relation. Overall, this approach reduces the computational complexity to O(NlogN). NOTE: Please include the following references in any work that utilizes this code: (1) Li, P., Johnston, H., Krasny, R.: A Cartesian Taylor series treecode for screened Coulomb interactions. J. Comput. Phys. 228 (2009) 3858-3868 (2) Duan, Z.-H., Krasny, R.: An adaptive treecode for computing nonbonded potential energy in classical molecular systems. J. Comput. Chem. 22 (2001) 184-195 (3) Lindsay, K., Krasny, R.: A particle method and adaptive treecode for vortex sheet motion in 3-D flow. J. Comput. Phys. 172 (2001) 879-907 (4) 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 ./source, ./example and ./aux : ----------------------------------------- ./source treedriver_yukawa.f : Driver program for testing the treecode subroutine treecode3d_yukawa.f. treecode3d_yukawa.f : Treecode subroutine. compile : Compile line for Intel ifort compiler. ./example xyz.dat : ASCII 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]. q.dat : ASCII file containing 50,000 charges uniformly distributed in [-0.5,0.5]. input_50k_p : Input data file provided for testing the treecode for computing the potential energy. The files xyz.dat and q.dat provide the particle data. output_50k_p : Program output using input_50k_p. These results were produced using the Intel ifort compiler and a single core of an Apple Mac Pro with dual 2.8GHz Harpertown quad-core CPUs and 16GB of memory. input_50k_pf : Input data file provided for testing the treecode for computing the potential energy AND force. The files xyz.dat and q.dat provide the particle data. output_50k_pf : Program output using input_50k_pf. These results were produced using the Intel ifort compiler and a single core of an Apple Mac Pro with dual 2.8GHz Harpertown quad-core CPUs and 16GB of memory. ./aux gendata3d.f : A program 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 xyz.dat and q.dat file. The user has the option of generating particles uniformly distributed in [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5], or these particle positions projected radially onto the surface of the unit sphere. Compiling and Testing the code ------------------------------ Compiling : The source code is written in standard Fortran90, and successfully compiled on a MAC OS X 10.5.7 with the Intel ifort compilers, as well as a RHEL 4.0 Linux machine with the Pathscale pathf90 compiler: Intel MAC : --------- ifort -fast treedriver_yukawa.f treecode3d_yukawa.f AMD Opteron Linux box : --------------------- pathf90 -Ofast treedriver_yukawa.f treecode3d_yukawa.f Testing : Once compiled, to test the code cd to the example directory and type ../source/a.out < ./input_50k_p This executes the code with N=50,000 and 4 other user specified parameters (see the description of inputs below). The relative error in both the L2 and infinity norms of the potential energy are output, as well as the total time for the calculation. One may also test the code using input_50k_pf in which iflag=2, in which case both the potential energy and force are computed. Description of inputs and treecode3d_yukawa.f calling statement: --------------------------------------------------------------------- treedriver_yukawa.f : ------------------- The program treedriver will prompt the user for the following: NUMPARS : (INTEGER) Number of particles KAPPA : (REAL*8) Debye-Huckel screening parameter 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_yukawa.f.) ORDER : (INTEGER) Order of 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. IFLAG : (INTEGER) Flag determining what is to be computed: IFLAG = 1 : potential energy IFLAG = 2 : potential energy AND forces treecode3d_yukawa.f : ------------------- SUBROUTINE TREECODE3D_YUKAWA(x,y,z,q,kappa,order,theta,maxparnode, & numpars,orderind,iflag,tpoten, & tforce,forcedim) CALL TREECODE3D_YUKAWA(x,y,z,q,kappa,order,theta,maxparnode,numpars, & orderind,iflag,tpoten,tforce,forcedim) Input: ----- X,Y,Z : (REAL*8 one-dimensional arrays of length NUMPARS) Cartesian coordinates of the particles. Q : (REAL one-dimensional array of length NUMPARS) Corresponding charge of each particle. KAPPA : (REAL*8) Debye-Huckel screening parameter 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. 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. NUMPARS : (INTEGER) number of particles IFLAG : (INTEGER) Flag determining what is to be computed: IFLAG = 1 : potential energy IFLAG = 2 : potential energy AND forces FORCEDIM : (INTEGER) leading dimension of TFORCE array - TFORCE(FORCEDIM,3) **NOTE** : If you only need the potential energy set IFLAG=1. You MUST still pass the TFORCE array. However, in this case it need only have one row, e.g. REAL(KIND=r8),DIMENSION(1,3) :: tforce in which case you should set FORCEDIM=1. Output: ------ TPOTEN : (REAL one-dimensional array of length NUMPARS) ith entry of the arrray is the approximate potential at the (reordered) ith particle TFORCE : (REAL two-dimensional array of size FORCEDIMx3) ith entry (i.e. (i,1:3)) of the arrray is the approximate force at the (reordered) ith particle ORDERIND : (INTEGER one-dimensional array of length NUMPARS) orderind(i) is new position/index of the particle that was initially in the ith position on entry to the subroutine. The reordering results from building the tree structure.