cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  This is the 3D fast multipole method (FMM) code for the Yukawa 
c  potential (also called modified Helmholtz, screened Coulomb, 
c  linearized Poisson-Boltzman, etc.). It is based on the new version 
c  of FMM first introduced by Greengard and Rokhlin for the Laplace
c  equation in 1997 (see Ref. 2 below). It computes the screened 
c  Coulombic interaction between n particles.
c
c  Copyright (C) Jingfang Huang 
c
c  This program is free software; you can redistribute it and/or
c  modify it under the terms of the GNU General Public License
c  as published by the Free Software Foundation; either version 2
c  of the License, or (at your option) any later version.
c
c  This program is distributed in the hope that it will be useful,
c  but WITHOUT ANY WARRANTY; without even the implied warranty of
c  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c  GNU General Public License for more details.
c
c  License details are available in license.txt, you may also write to 
c
c    The Free Software Foundation, Inc., 
c    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
c
c  Reference:
c    1. Greengard, Leslie F.; Huang, Jingfang: A new version of the 
c       fast multipole method for screened Coulomb interactions in 
c       three dimensions. J. Comput. Phys. 180 (2002), no.2, 642--658. 
c 
c    2. Greengard, Leslie; Rokhlin, Vladimir: A new version of the 
c       fast multipole method for the Laplace equation in three dimensions.  
c       Acta numerica, 1997,  229--269, Acta Numer., 6, Cambridge Univ. 
c       Press, Cambridge, 1997. 
c
c  For suggestions, comments, and bug reports, please contact 
c
c    Jingfang Huang
c    CB# 3250, Phillips Hall
c    Department of Mathematics, UNC
c    Chapel Hill, NC 27599-3250, USA.
c    Email: huang@amath.unc.edu
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  this is the driver of the uniform fast multipole code.
c
c  the following variables should be defined in the driver.
c
c    natoms: total number of atoms.
c    beta: the frequency of the equation.
c    zat(3,natoms): particle locations.
c    charge(natoms) : charge each particle carries.
c    pot(natoms): potential at each particle location.
c    field(3,natoms): field at each particle location.
c    nlev: number of levels in the octree structure.
c    ier: error message.
c
c  a subroutine is provided to compare fmm results with direct results,
c  so you also provide the following vectors for such comparisons.
c
c    dpot(natoms): potential using direct summation.
c    dfield(3,natoms): field using direct summation.
c
c  for advanced users:
c    if you want to change the precision of this code, you can check
c    parm-uniyuk.h.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NATOMS,NLEV,IER
      PARAMETER (NATOMS=100000)
      REAL *8 BETA
      REAL *8 ZAT(3,NATOMS),CHARGE(NATOMS)
      REAL *8 POT(NATOMS),DPOT(NATOMS)
      REAL *8 FIELD(3,NATOMS),DFIELD(3,NATOMS)
c
c-----functions called
c
      REAL *8 SECOND
c
c-----local variables, mostly for comparing direct results with fmm
c     results.
c
      INTEGER *4 I,IMAX
      REAL *8 CTOT,TIME0,TIME1,SALG,STOT,SALG2,STOT2
      REAL *8 TT,EALG,ERRMAX
c
c-----call initialization routines
c     1. initialize the output units.
c        6: screen.
c        13: fort.13.
c
      CALL PRINI(6,13)
c
c-----2. set up parameters.
c
      BETA=0.1D0
      NLEV=7
c
      CALL PRIN2('frequency = *', BETA, 1)
      CALL PRINF('nlev = *',NLEV,1)
      CALL PRINF('natoms = *',NATOMS,1)
c
c-----3. generate charges and their locations.
c        for uniform code, all particles should be located
c        insdie the unit box. i.e., for other distributions
c        you have to rescale your box size.
c
      CALL DUMMY(NATOMS,ZAT,CHARGE)
c
c======================================================================d
c-----4. main fmm call, call fmm to calculate the potential
c
      CALL FMMYUK_UNI(BETA,NATOMS,ZAT,CHARGE,POT,FIELD,NLEV,IER)
c
c-----end of the fast multipole algorithm.
c======================================================================d
c
c-----5. finally, comparison with direct method.
c       write out first imax data points
c
      IMAX = NATOMS
c      IF (NATOMS .GT. 400) IMAX = 400
c
c      call prin2( 'from alg, pot = *',pot,imax)
c      call prin2( 'from alg, field = *',field,3*imax)
c
c-----the direct method. parameter imax below determines how many
c     direct interactions are computed.
c
      TIME0 = SECOND()
      DO 1100 I=1,NATOMS
        CALL YDIRECI(BETA,NATOMS,ZAT,CHARGE,I,DPOT(I),DFIELD(1,I))
1100  CONTINUE
      TIME1 = SECOND()

      OPEN (UNIT = 18 , FILE = 'points_in_space2', FORM = 'unformatted',
     1      access = 'direct' , recl = 8*3*100000)
      OPEN (UNIT = 19 , FILE = 'direct_potential', FORM = 'unformatted',
     1      access = 'direct' , recl = 8*100000) 
      OPEN (UNIT = 20 , FILE = 'direct_field', FORM = 'unformatted',
     1      access = 'direct' , recl = 8*3*100000) 
c      OPEN (UNIT = 21 , FILE = 'direct_test_output')


      WRITE( 18 , rec = 1 ) ZAT
      WRITE( 19 , rec = 1 ) DPOT
      WRITE( 20 , rec = 1 ) DFIELD
      CLOSE (20)
      CLOSE (19)
      CLOSE (18)
     
      OPEN (UNIT = 21 , FILE = 'direct_test_output')
      DO 2500 I = 1, NATOMS            
      WRITE (21,*) ZAT(1,I) ,ZAT(2,I) ,ZAT(3,I) ,DPOT(I) ,
     1      DFIELD(1,I) ,DFIELD(2,I) ,DFIELD(3,I) 
2500  CONTINUE

      CLOSE(21)


c
c-----error analysis.
c
c      SALG = 0.0D0
c      SALG2= 0.0D0
c      STOT = 0.0D0
c      STOT2= 0.0D0
c      ERRMAX=0.0D0
c      DO 1200 I=1,IMAX
c        SALG = SALG + (DPOT(I) - POT(I))**2
c        STOT = STOT + DPOT(I)*DPOT(I)
c        SALG2= SALG2+ (DFIELD(1,I)-FIELD(1,I))**2+
c     1    (DFIELD(2,I)-FIELD(2,I))**2+(DFIELD(3,I)-FIELD(3,I))**2
c        STOT2= STOT2+ DFIELD(1,I)*DFIELD(1,I)+
c     1    DFIELD(2,I)*DFIELD(2,I)+DFIELD(3,I)*DFIELD(3,I)
c        IF (DABS(DPOT(I)-POT(I)).GT.ERRMAX) THEN
c          ERRMAX=DABS(DPOT(I)-POT(I))
c        ENDIF
1200  CONTINUE
c
c      CALL PRIN2(' time for imax points directly is *',TIME1-TIME0,1)
c      TT = (TIME1-TIME0)*NATOMS/IMAX
      TT = (TIME1-TIME0)
      CALL PRIN2(' time for all points directly is *',TT,1)

      WRITE(11,556)TT
556   FORMAT(' time for direct calc is ',F12.3)
c      call prin2( 'direct pot = *',dpot,imax)
c      call prin2( 'direct field = *',dfield,3*imax)
c      EALG = DSQRT(SALG/STOT)
c      CALL PRIN2( ' l2 error of fast alg = *',EALG,1)
c      CALL PRIN2( ' max error of potential = *',ERRMAX,1)
c      EALG = DSQRT(SALG2/STOT2)
c      CALL PRIN2( ' l2 error of field = *',EALG,1)
c
      STOP
      END
c
