Masayuki Toyoda and Taisuke Ozaki

*Research Center for Integrated Science,
Japan Advanced Institute of Science and Technology,
*

1-1, Asahidai, Nomi, Ishikawa 932-1292, Japan

1-1, Asahidai, Nomi, Ishikawa 932-1292, Japan

NOTICE: This is the authorfs version of a work accepted for publication by Elsevier. Changes resulting from the publishing process, including peer review, editing, corrections, structural formatting and other quality control mechanisms, may not be reflected in this document. Changes may have been made to this work since it was submitted for publication. A definitive version was subsequently published in Computer Physics Communications

Abstract

We propose a new method for the numerical evaluation of the spherical Bessel transform. A formula is derived for the transform by using an integral representation of the spherical Bessel function and by changing the integration variable. The resultant algorithm consists of a set of the Fourier transforms and numerical integrations over a linearly spaced grid of variable in Fourier space. Because the -dependence appears in the upper limit of the integration range, the integrations can be performed effectively in a recurrence formula. Several types of atomic orbital functions are transformed with the proposed method to illustrate its accuracy and efficiency, demonstrating its applicability for transforms of general order with high accuracy.

Keywords: Hankel transforms, spherical Bessel functions, atomic orbitals

PACS: 02.30.Uu, 71.15.Ap, 31.15.-p

PACS: 02.30.Uu, 71.15.Ap, 31.15.-p

Introduction

Several computational methods have been developed for the Hankel transforms [6,7,8]. However, not many of them can be applied for SBT. There have been proposed three different approaches for SBT for general order: (a) recasting the transform integral as a convolution integral by changing the coordinate variables [9,10,11], (b) expansion in terms of a series of Fourier cosine and sine transforms by the trigonometric expansion of the spherical Bessel function [12], and (c) the discrete Bessel transform method which describes SBT as an orthogonal transform [13]. The approach (a) is quite fast since it utilizes the fast Fourier transform (FFT) algorithm. The computation time actually scales as , where is the number of quadrature points. It, however, has a major drawback caused by the logarithmic grid that almost all the grid points are located in close proximity of the origin. The computation time of the approach (b) which also uses FFT scales as , where is the order of transform. Therefore, it becomes slower for the higher order transform. In addition to that, since it requires the integrand to be multiplied by inverse powers of the radial coordinate, the high order transforms may become unstable. The computation time of the approach (c) scales as . The quadrature points are located at each zero of the spherical Bessel function. The optimized selection of the quadrature points enables us to use a small number of while keeping the accuracy of the computation. However, when consecutive transforms with different orders are required, it may become a minor trouble that the optimized quadrature points differ depending on the order of transform.

In this paper, we propose a new method for the numerical SBT which uses a linear coordinate grid. The transform is decomposed into the Fourier transforms and the numerical integrations which can be evaluated recursively. The computation time for the present method scales as with overhead for the numerical integration which scales as . The linear coordinate grid prevents us from troubles caused by the non-uniform or order-dependent grid points. If the considered problem requires to transform a function with various orders, the present method has further the advantage that the results of the most time consuming calculations (i.e. the Fourier transforms and the integrations) for a transform with a certain order can also be used for transforms with any order less than .

Formulation

where is the spherical Bessel function of the first kind. The integral representation of the spherical Bessel function is given by

where is the Legendre polynomials

(3) | ||

(4) |

Here, is the double factorial

(5) |

with an exceptional definition that . By substituting Eq. (2) into Eq. (1), and using the parity property , the transform is rewritten as follows:

(6) |

Now, we change the variables as , and it becomes

(7) |

where is the Fourier cosine/sine transform of ,

and is the largest integer that does not exceed . Finally, by expanding the Legendre polynomials, the transform is decomposed into a sum of definite integrals

where is given by

Since the integrals appearing in Eqs. (9) and (10) have the same parity as the order of transform , either the Fourier cosine or sine transform is required to be performed, for given . In our implementation, as explained later, at most two more Fourier transforms are required to be performed to evaluate the derivatives of . Therefore, regardless of the order of transform, only three Fourier transforms are required. On the other hand, since the integrals depend on the order of transform through terms, a number of different integrals are required for the summation in Eqs. (9) and (10). Even so, however, this does not increase the computational cost because does not appear in the integrand of Eq. (11) and thus the computation cost for scales as , rather than as .

Implementation

(12) | ||

(13) |

At each point of , the integral is divided into its segments

where the segments are defined as

(15) | ||

(16) |

By defining the sum of the segments as , the following simple recurrence formula is obtained:

Therefore, the evaluation of the integral is accomplished for all the -grid points through a summation of segments , where, at each step of the summation, the subtotal divided by gives .

Each segment is evaluated by locally interpolating with a polynomial curve. Since no grid point is available in between the both ends of the integration range, the derivatives of are also required to interpolate with higher order polynomials. By noting that only the trigonometric functions depend on in the integrand of (8), the derivatives and second derivatives are obtained analytically as follows:

By interpolating with a -th order polynomial, the integral segments are given by

where the coefficients are determined from the values and derivatives of at and (see Appendix A). We have performed the interpolation with linear (), cubic (), and quintic () polynomials. Since a higher order interpolation suffers more severely from the arithmetic errors, interpolations with a polynomial of order have not been tested.

To summarize, the present method computes the transform of a function with the order via the following steps:

- Fourier cosine/sine transforms of multiplied by a power of (Eqs. (8), (19), and (20)).
- Evaluation of the integral segments by the piecewise polynomial interpolation. (Eq. (21)).
- Evaluation of through the summation of (Eqs. (17) and (18)).
- Weighted summation of to give the transformed function (Eqs. (9) and (10)).

Computation results

The normalized GTO function in spherical coordinate is defined as follows:

(22) |

where the normalization factor is given by

(23) |

The corresponding transformed function is

(24) |

Figure 1 (a) shows the the numerical error in the present method in the 0th-order transform of the GTO function, where the integral segments have been evaluated by interpolating with the linear (solid line), cubic (dashed line), and quintic (dash-dotted line) polynomials. The parameters used in the calculation are as follows: the exponent of GTO , the number of grid points , and the maximum value of -grid . It is clearly observed that the error is reduced quickly as the order of interpolation polynomial increases, and that the quintic polynomial gives a sufficiently accurate result. In Fig. 1 (b), the numerical error for the 15th-order transform of GTO is shown, where the same parameters as the previous calculation and the quintic polynomial interpolation is used. Even in such a high order transform, the error remains small. This implies that the numerical error comes mainly from the integration of the segments while the possible rounding-off error in Eqs. (9), (10), and (11) is actually negligible unlike in the approach (b) referred in the introduction of this paper [12].

Since the numerical error in the integration of a segment remains in the summation in Eq. (14), the numerical error of the segments in the small- region also contributes to the error in the large- region. This explains why the damping of the error at large- region is so slow in Fig. 1 (a). A downward summation is effective to reduce this error. In Fig. 2, the numerical error in the 0th-order transform of GTO is plotted, where the summation of the segments are performed upward (solid line) and downward (dashed line). In the downward summation, the numerical error in the large- region becomes much smaller, while, in the small- region, the error becomes larger because of the accumulation of the error from the large- region. Therefore, by connecting the results of the upward and downward summations at a certain point (for example, = 5), accurate results can be obtained in both small- and large- regions.

So far, the transforms have been performed where the order of transform and the order of the GTO function are equivalent. In Fig. 3, the numerical error in the transform with the order for the GTO function whose order is is plotted with a variety of grid spacings in real space. It is found that the error (solid line) is larger than that of the case calculated with the same condition (dash-dotted line in Fig. 1 (a)). The decrease of the error by smaller grid spacing implies that fine grid spacing is necessary for the accurate sine transform in Eq. (8) since the GTO function of the order has a finite value at while the spherical Bessel function of the order vanishes proportionally to . The different behaviors near the origin results in another undesirable fact that the transformed function in Fourier space has a very long tail. In the case considered here, the analytic form of the SBT of the GTO function is given as

(25) |

where is the Dawson's integral [14] which is defined as

(26) |

From the asymptotic form of the Dawson's integral, the SBT of the GTO function is known to behave as

(27) |

when is large. Therefore, in general, a wide range of -grid has to be employed where the orders of transform and the function are not equivalent.

In quantum chemistry, the STO wave functions are often used as the basis functions for the molecular orbitals. The radial part of the STO functions is given by

(28) |

where the normalization factor is

(29) |

The corresponding transformed function is given as

(30) |

In Fig. 4, plotted is the numerical error for the 0th- and 15th-order transforms of STO. The parameters used in the calculations are as follows: the number of grid points and the maximum value of the -grid (solid line in Fig. 4 (a)); and (dashed line in Fig. 4 (a)); and (solid line in Fig. 4 (b)); and (dashed line in Fig. 4 (b)). The exponent of STO is for all the calculations. There are two kinds of sources for the numerical error in the SBT for STO: the cusp at the origin and the long tail of STO, and they are illustrated in Fig. 4. The use of the fine grid by increasing reduces the error as shown in Fig. 4 (a), which implies the requirement of a fine grid near the cusp for the accurate integration in Eq. (8). On the other hand, the error is reduced by using the wide range of -grid as shown in Fig. 4 (b) because of the long tail of STO ().

In the electronic structure calculations based on the density functional theory (DFT), the non-analytic atomic orbitals are often used as well as the analytic functions such as the GTO and STO functions. The PAO function is one of those non-analytic basis functions and used in many DFT calculation methods. A PAO function is calculated by solving the atomic Kohn-Sham equation with confinement pseudopotentials [15]. Figure 5 (a) shows the PAO functions of the and states of an oxygen atom, where the confinement radius is a.u. The corresponding transformed functions are plotted in Fig. 5 (b). Since the mathematical formula for the forward and backward SBT is equivalent except for an additional factor of , the function should recover the original shape by two consecutive transforms and thus we can illustrate the numerical error by comparing the input function and the back-and-forth transformed function. In Fig. 6, the numerical error accompanied by the present method in two consecutive transforms of the PAO functions of the oxygen and states is plotted. The parameters used for the calculations are: and . The error is less than which may not be sufficiently small but is acceptable in practice depending on the particular problems.

For the purpose of comparison, we have performed the back-and-forth transform of the PAO function of the oxygen state with the method proposed by Siegman and Talman [9,10]. In the calculations, we have used and . As shown in Fig. 7, to suppress the error less than in the region , a very large number of quadrature points () are required. As in this case, due to the use of the linear coordinate grid, a small number of the grid points is enough for the present method to achieve the same degree of accuracy. The computation time has also been measured by taking the average of the cpu time used for 1000 repetitions of a set of transforms, which was carried out on a machine with an Intel Core i7 processor at 2.67 GHz. The program code of the Siegman-Talman method is also implemented by the authors according to the article [10]. The Siegman-Talman method costs 0.64 msec per transform with , and the present method costs 1.32 msec per transform with . In this measurement, each transform is performed independently so that any speed-up technique is not used such as skipping redundant calculations and reusing recources. This is just a rough estimation, but looks reasonable considering the theoretical computation cost of the both methods (i.e. the computation cost for the Siegman-Talman method is basically coming from that of two Fourier transforms, whereas the present method requires three times of the Fourier transforms as well as the additional numerical integrations). In the present method, as mentioned before, once a transform with an order is performed, another transform with another order can also be performed by repeating only the final step. In the set of transform with 15 different orders for a single input function, a substantial amount of redundant calculations are included. In fact, similar redundancy arises also in the Siegman-Talman method. By skipping such redundant calculations, average time per transform reduces to 0.45 msec for the Siegman-Talman method and to 0.22 msec for the present method, in the same measurement as previous one except for the skipping. Therefore, the present method can be quite effective for a particular type of applications where a single input function is required to be transformed several times with different orders. An example of such applications is the calculation of two-electron integrals in molecules and solids [4,5]. The authors believe that the present method is an alternative to the Siegman-Talman method with comparable efficiency. The most significant difference of the present method from the Siegman-Talman method is the use of a uniform coordinate grid. For many physical systems which are better described by the uniform coordinate grid, a small number of the quadrature points are required, which gives the advantages not only in the computation speed but also in, for example, the reduction in the memory size and communication traffic in parallelized computations.

Summary

Acknowledgments

(A.1) |

where the coefficients are

(A.2) | ||

(A.3) | ||

(A.4) | ||

(A.5) |

A quintic () polynomial is obtained by using the second derivatives of ,

(A.6) |

where the coefficients are

(A.7) | ||

(A.8) | ||

(A.9) | ||

(A.10) | ||

(A.11) | ||

(A.12) |

(B.1) | ||

(B.2) | ||

(B.3) | ||

(B.4) |

where

(B.5) |

Therefore, the segment for an even number of is given as

(B.6) | ||

(B.7) |

Similarly, for an odd number of , it is

(B.8) |

- [1]
L. W. Person, P. Benioff, Calculations for quasi-elastic scattering
,
and
at 460
, Nucl. Phys.
**A187**(1972) 401. - [2]
A. E. Glassgold, G. Ialongo, Angular distribution of the outgoing electrons in
electronic ionization, Phys. Rev.
**175**(1968) 151. - [3]
M. Liguori, A. Yadav, F. K. Hansen, E. Komatsu, S. Matarrese, B. Wandelt,
Temperature and polarization CMB maps from primordial non-Gaussiantities of
the local type, Phys. Rev. D
**76**(2007) 105016. - [4]
J. D. Talman, Numerical methods for multicenter integrals for numerically
defined basis functions applied in molecular calculations, Int. J. Quantum
Chem.
**93**(2003) 72. - [5]
M. Toyoda, T. Ozaki, Numerical evaluation of electron repulsion integrals for
pseudoatomic orbitals and their derivatives, J. Chem. Phys.
**130**(2009) 124114. - [6]
S. M. Candel, An algorithms for the Fourier Bessel transform, Comput. Phys. Comm.
**23**(1981) 343. - [7]
J. D. Secada, Numerical evaluation of the Hankel transform, Comput. Phys. Comm.
**116**(1999) 278. - [8]
V. K. Singh, Om P. Singh, R. K. Pandey, Numerical evaluation of the Hankel
transform by using linear Legendre multi-wavelets, Comput. Phys. Comm.
**179**(2008) 424. - [9]
A. E. Siegman, Quasi fast Hankel transform, Opt. Lett.
**1**(1977) 13. - [10]
J. D. Talman, Numerical Fourier and Bessel transforms in logarithmic variables,
J. Comput. Phys.
**29**(1978) 35. -
J. D. Talman, NumSBT: A subroutine for calculating spherical Bessel transforms
numericaly, Computer Phys. Comm.
**180**(2009) 332. - [12]
O. A. Sharafeddin, H. F. Bowen, D. J. Kouri, D. K. Hoffman, Numerical
evaluation of spherical Bessel transforms via fast Fourier transforms, J.
Comput. Phys.
**100**(1992) 294. - [13]
D. Lemoine, The discrete Bessel transform algorithm, J. Chem. Phys.
**101**(1994) 3936. - [14]
M. Abramowitz, I. A. Stegun (Eds.),
*Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables*, Dover publications, 1972. - [15]
T. Ozaki, Variationally optimized atomic orbitals for large-scale electronic
structures, Phys. Rev. B
**67**(2003) 155108.

This document was generated using the
**LaTeX**2`HTML` translator Version 2002-2-1 (1.71)

Copyright © 1993, 1994, 1995, 1996,
Nikos Drakos,
Computer Based Learning Unit, University of Leeds.

Copyright © 1997, 1998, 1999,
Ross Moore,
Mathematics Department, Macquarie University, Sydney.

TOYODA Masayuki 2009-09-30