Error estimates for the fast multipole method

Error estimates for algorithms based on truncations for evaluating electrostatic interactions in molecular dynamics applications are very important for several reasons. For example, the estimates are necessary to establish the validity of the simulations and can be used to estimate various simulation parameters. Very precise estimates have been found for the Ewald method and the related particle mesh Ewald method. However, for the very popular fast multipole method such a precise estimate is not available. In this paper, we illustrate the rather complicated error behavior of the fast multipole method and we use statistical methods to derive an estimate for the root mean square error on the forces. Furthermore, the expected maximum error on the force acting on a single particle is studied. The estimates are tested against errors obtained from simulations and are found to be very precise.     

Implementing the fast multipole method in three dimensions

The Rokhlin-Greengard fast multipole algorithm for evaluating Coulomb and multipole potentials has been implemented and analyzed in three dimensions. The implementation is presented for bounded charged systems and systems with periodic boundary conditions. The results include timings and error characterizations.     

多层快速多极子算法中的两步插值技术

多层快速多极子算法(MLFMA)在快速多极子算法(FMM)的基础上按多层聚集、层间转移和多层扩散的思路以达到优化矩阵向量积的运算的目的,其中多层聚集和多层扩散过程,随着层数递增,角谱积分采样点数逐层递增,为了快速计算角谱积分,需要采用插值技术和反插值技术以提高计算效率.应用两步插值技术替代传统的单步插值技术,大幅提高了多层快速多极子层间插值反插值操作的计算效率,对于应用普通个人计算机求解特大电大尺寸问题,具有重要意义.     

Acceleration of canonical molecular dynamics simulations using macroscopic expansion of the fast multipole method combined with the multiple timestep integrator algorithm

A canonical molecular dynamics (MD) simulation was accelerated by using an efficient implementation of the multiple timestep integrator algorithm combined with the periodic fast multiple method (MEFMM) for both Coulombic and van der Waals interactions. Although a significant reduction in computational cost has been obtained previously by using the integrated method, in which the MEFMM was used only to calculate Coulombic interactions (Kawata, M., and Mikami, M., 2000, 98, J. Comput. Chem., in press), the extension of this method to include van der Waals interactions yielded further acceleration of the overall MD calculation by a factor of about two. Compared with conventional methods, such as the velocity-Verlet algorithm combined with the Ewald method (timestep of 0.25 fs), the speedup by using the extended integrated method amounted to a factor of 500 for a 100 ps simulation. Therefore, the extended method reduces substantially the computational effort of large scale MD simulations.     

Spectral accuracy in fast Ewald-based methods for particle simulations

A spectrally accurate fast method for electrostatic calculations under periodic boundary conditions is presented. We follow the established framework of FFT-based Ewald summation, but obtain a method with an important decoupling of errors: it is shown, for the proposed method, that the error due to frequency domain truncation can be separated from the approximation error added by the fast method. This has the significance that the truncation of the underlying Ewald sum prescribes the size of the grid used in the FFT-based fast method, which clearly is the minimal grid. Both errors are of exponential-squared order, and the latter can be controlled independently of the grid size. We compare numerically to the established SPME method by Essmann et al. and see that the memory required can be reduced by orders of magnitude. We also benchmark efficiency (i.e. error as a function of computing time) against the SPME method, which indicates that our method is competitive. Analytical error estimates are proven and used to select parameters with a great degree of reliability and ease.     

边界元奇异与近奇异数值积分方法及其应用于大规模声学问题

提出了综合处理Burton-Miller方法所导致的奇异积分与近奇异积分问题的数值求积方法,以此改进了基于常量元素的常规边界元和低频快速多极边界元方法。对于奇异积分问题,利用Hadamard有限积分方法进行解决;对于近奇异积分问题,则采用极坐标变换法和PART方法(Projection and Angular&;Radial Transformation)进行克服。与解析解和LMS Virtual.Lab商业软件的结果比较验证了方法的正确性,并对比分析了奇异积分与近奇异积分对计算精度的影响。采用低频快速多极子方法以加速常规边界元法的计算效率,计算分析了计算复杂度,并成功实现了34万自由度大规模问题的计算。结果表明,近奇异积分问题主要由超奇异核函数引起,对计算精度的影响不容忽略;快速多极边界元法的精度与常规边界元法一致,但计算复杂度要远低于后者。      

Numerical quadrature for singular and near-singular integrals of boundary element method and its applications in large-scale acoustic problems

The numerical quadrature methods for dealing with the problems of singular and near-singular integrals caused by Burton-Miller method are proposed,by which the conventional and fast multipole BEMs(boundary element methods) for 3D acoustic problems based on constant elements are improved.To solve the problem of singular integrals,a Hadamard finite-part integral method is presented,which is a simplified combination of the methods proposed by Kirkup and Wolf.The problem of near-singular integrals is overcome by the simple method of polar transformation and the more complex method of PART(Projection and Angular Radial Transformation).The effectiveness of these methods for solving the singular and near-singular problems is validated through comparing with the results computed by the analytical method and/or the commercial software LMS Virtual.Lab.In addition,the influence of the near-singular integral problem on the computational precisions is analyzed by computing the errors relative to the exact solution.The computational complexities of the conventional and fast multipole BEM are analyzed and compared through numerical computations.A large-scale acoustic scattering problem,whose degree of freedoms is about 340,000,is implemented successfully.The results show that,the near singularity is primarily introduced by the hyper-singular kernel,and has great influences on the precision of the solution.The precision of fast multipole BEM is the same as conventional BEM,but the computational complexities are much lower.     

新的对角形式快速多极边界元法求解声学Helmholtz方程

传统外部声学Helmholtz边界积分方程无法在个人计算机上求解大规模工程问题. 为了有效解决这个问题, 将快速多极方法引入到边界积分方程中, 加速系统矩阵方程组的迭代求解. 由于在边界积分方程中引入基本解的对角形式多极扩展, 新的快速多极边界元法的计算效率与传统边界元相比显著提高, 计算量和存储量减少到O(N)量级(N为问题的自由度数). 包括含有420000个自由度的大型潜艇模型数值算例验证了快速多极边界元法的准确性和高效性, 清楚表明新算法在求解大规模声学问题中的优势,     

新的对角形式快速多极边界元法求解声学Helmholtz方程

传统外部声学Helmholtz边界积分方程无法在个人计算机上求解大规模工程问题. 为了有效解决这个问题, 将快速多极方法引入到边界积分方程中, 加速系统矩阵方程组的迭代求解. 由于在边界积分方程中引入基本解的对角形式多极扩展, 新的快速多极边界元法的计算效率与传统边界元相比显著提高, 计算量和存储量减少到O(N)量级(N为问题的自由度数). 包括含有420000个自由度的大型潜艇模型数值算例验证了快速多极边界元法的准确性和高效性, 清楚表明新算法在求解大规模声学问题中的优势, 具有良好的工程应用前景.     

On solving continuum-mechanics problems by fast multipole methods

It is shown analytically and numerically that the kernel-independent fast multipole method provides the accuracy and computational complexity of the analytical method if circular (spherical in 3D) equivalent surfaces are used.     

Computation of scattering from N spheres using multipole reexpansion

A computational technique for the solution of problems of wave scattering from multiple spheres is developed. This technique, based on the T-matrix method, uses the theory for the translation and reexpansion of multipole solutions of the Helmholtz equation for fast and exact recursive computation of the matrix elements. The spheres can have prescribed radii, impedances, and locations. Results are validated by comparison with boundary element calculations, and by convergence analyses. The method is much faster than numerical methods based on discretization of space, or of the sphere surfaces. An even faster method is presented for the case when the spheres are aligned coaxially.     

快速多极边界元法用于扩散光学断层成像研究

在求解扩散光学断层成像中的正向问题时, 目前普遍采用有限元法, 但是随着实际模型规模的增大, 有限元法的计算量问题日益显著, 而边界元法则由于可以降低计算维度使计算量减少而备受关注. 本文以均匀的高散射介质为模型, 研究了将快速多极边界元法用于扩散光学断层成像的正向问题. 快速多极边界元法利用核函数的多极展开, 将常规边界元法中系数矩阵和迭代矢量的乘积项等价为相应四叉树结构的一次递归, 再结合广义最小残量法进行迭代求解. 将计算结果和蒙特卡罗法的模拟结果进行了比较, 表明利用快速多极边界元法的模拟结果和蒙特卡罗法的结果有很好的一致性. 研究结果验证了快速多极边界元法可以用于扩散光学断层成像, 为其大规模和实时成像带来可观的前景. 关键词： 扩散光学断层成像 边界元法 快速多极边界元法     

Convergence analysis of the Wolf method for Coulombic interactions

A rigorous proof for convergence of the Wolf method (Wolf et al., 1999 [9]) for calculating electrostatic energy of a periodic lattice is presented. In particular, we show that for an arbitrary lattice of unit cells, the lattice sum obtained via Wolf method converges to the one obtained via Ewald method.     

Preconditioning based on Calderon’s formulae for periodic fast multipole methods for Helmholtz’ equation

Solution of periodic boundary value problems is of interest in various branches of science and engineering such as optics, electromagnetics and mechanics. In our previous studies we have developed a periodic fast multipole method (FMM) as a fast solver of wave problems in periodic domains. It has been found, however, that the convergence of the iterative solvers for linear equations slows down when the solutions show anomalies related to the periodicity of the problems. In this paper, we propose preconditioning schemes based on Calderon’s formulae to accelerate convergence of iterative solvers in the periodic FMM for Helmholtz’ equations. The proposed preconditioners can be implemented more easily than conventional ones. We present several numerical examples to test the performance of the proposed preconditioners. We show that the effectiveness of these preconditioners is definite even near anomalies.     

A fast Fourier transform on multipoles (FFTM) algorithm for solving Helmholtz equation in acoustics analysis

This article presents a fast algorithm for the efficient solution of the Helmholtz equation. The method is based on the translation theory of the multipole expansions. Here, the speedup comes from the convolution nature of the translation operators, which can be evaluated rapidly using fast Fourier transform algorithms. Also, the computations of the translation operators are accelerated by using the recursive formulas developed recently by Gumerov and Duraiswami [SIAM J. Sci. Comput. 25, 1344-1381(2003)]. It is demonstrated that the algorithm can produce good accuracy with a relatively low order of expansion. Efficiency analyses of the algorithm reveal that it has computational complexities of O(Na), where a ranges from 1.05 to 1.24. However, this method requires substantially more memory to store the translation operators as compared to the fast multipole method. Hence, despite its simplicity in implementation, this memory requirement issue may limit the application of this algorithm to solving very large-scale problems.     

二维声场预测的快速多极基本解法

传统基本解法在二维大规模模型的声场求解过程中,系统方程形成和求解的计算量正比于自由度N的二次方O(N2)和三次方O(N3),求解效率低;为此,引入快速多极子算法并采用广义极小残差法迭代求解,提出一种用于二维声场预测的快速多极基本解法。对无限长圆柱体及二维类车体辐射模型的仿真结果表明,当N为3000时,分别采用快速多极基本解法与传统基本解法求解所需的时间比值约为百分之四,且N越大比值越小;最终实现系统方程的形成和求解的计算量降低到正比于自由度O(N),提高了对二维大规模模型声场预测计算效率。      

大规模宽频弹性动力学分析的快速定向压缩边界元法

发展一种大规模宽频弹性动力学分析的快速定向压缩边界元法.证明弹性动力学核函数具有定向低秩特性，为采用快速定向压缩算法提供理论基础.根据S波的波数，将节点之间的相互作用划分为低频相互作用和高频相互作用，并将后者进一步划分为与多个楔形区.在楔形区上，可以采用核函数的定向低秩特性进行快速计算.低频相互作用与核无关快速多极边界元法中计算方法相同，不同方向楔形区上的变换矩阵可以采用坐标系旋转的方法进行快速计算.可对任意频率进行快速谐响应分析.数值算例表明：该方法可以将宽频弹性动力学问题计算复杂度降低到O（N log^αN）.与卷积求积边界元法相结合，也可以应用于弹性动力学瞬态分析.     

Yukawa potentials in systems with partial periodic boundary conditions. I. Ewald sums for quasi-two-dimensional systems

Yukawa potentials are often used as effective potentials for systems such as colloids, plasmas, etc. When the Debye screening length is large, the Yukawa potential tends to the non-screened Coulomb potential; in this small screening limit, or Coulomb limit, the potential is long-ranged. As is well known in computer simulation, a simple truncation of the long-ranged potential and the minimum image convention are insufficient to obtain accurate numerical data on systems. The Ewald method for bulk systems, i.e. with periodic boundary conditions in all three directions of space, has already been derived for the Yukawa potential [Molec. Phys. 88, 1357 (1996); J. Chem. Phys. 113, 10459 (2000)], but for systems with partial periodic boundary conditions, the Ewald sums have only recently been obtained [J. Chem. Phys. 126, 056101 (2007)]. In this paper, we provide a closed derivation of the Ewald sums for Yukawa potentials in systems with periodic boundary conditions in only two directions and for any value of the Debye length. Special attention is paid to the Coulomb limit and its relation to the electroneutrality of systems.     

Modeling and optimization of non-periodic grating couplers

The higher degree of freedom available for non-periodic gratings (as compared with their periodic counterparts) is investigated. These non-periodic structures may be employed to design novel light couplers with increased functionality. Optimizing such devices requires a complex search in a huge parameter space. The success in the solution of this task depends on the availability of a fast forward solver and a reliable search algorithm. Here, a fast forward solver based on the multiple multipole (MMP) method together with a near-to-far field transformation and a multiple scattering calculation is presented. Thanks to the efficiency of our approach, non-periodic gratings are evaluated with a speed comparable to commonly used periodic grating approximations. This allows our solver to be combined with a heuristic global search scheme, namely an evolutionary algorithm. The procedure is demonstrated with the optimization of a non-periodic grating output coupler that suppresses an unwanted second diffracted order.     

On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions

Methods for the evaluation of the electrostatic potential in systems which are periodic in three dimensions, where the electronic charge distribution is expanded as a linear combination of Gaussian functions of arbitrary quantum number, are described. An ‘exact’ method based on the Ewald potential function is derived, and an approximate scheme using a distributed point multipole representation of the charge distribution is then presented. Recursion formulae enable Cartesian and thence spherical derivatives (of any order) of the potential to be computed. Related procedures for the evaluation of the Fock matrix elements and the coulombic contribution to the total energy per unit cell are described. Tests of the exact and approximate procedures establish their relative accuracy and cost in the MgO, Si and Si₁₂O₂₄ systems in the rock-salt, diamond and chabazite structures, respectively. Applications include a study of the electric field gradient at the C, N and O sites in urea (crystalline and molecular), and contour mapping of the electrostatic potential in MgO and Si₁₂O₂₄.