Search

In this article we present an algorithm to efficiently evaluate the exchange matrix in periodic systems when Gaussian basis set with pseudopotentials are used. The usual algorithm for evaluating exchange matrix scales cubically with the system size because one has to perform O(N2) fast Fourier transforms (FFT). Here we introduce an algorithm that retains the cubic scaling but reduces the prefactor significantly by eliminating the need to do FFTs during each exchange build. This is accomplished by representing the products of Gaussian basis function using a linear combination of an auxiliary basis the number of which scales linearly with the size of the system. We store the potential due to these auxiliary functions in memory which allows us to obtain the exchange matrix without the need to do FFT, albeit at the cost of additional memory requirement. Although the basic idea of using auxiliary functions is not new, our algorithm is cheaper due to a combination of three ingredients: (a) we use robust Pseudospectral method that allows us to use a relatively small number of auxiliary basis to obtain high accuracy (b) we use occ-RI exchange which eliminates the need to construct the full exchange matrix and (c) we use the (interpolative separable density fitting) ISDF algorithm to construct these auxiliary basis that are used in the robust pseudospectral method. The resulting algorithm is accurate and we note that the error in the final energy decreases exponentially rapidly with the number of auxiliary functions.

We consider fast deterministic algorithms to identify the “best” linearly independent terms in multivariate mixtures and use them to compute, up to a user-selected accuracy, an equivalent representation with fewer terms. Importantly, the multivariate mixtures do not have to be a separated representation of a function. One algorithm employs a pivoted Cholesky decomposition of the Gram matrix constructed from the terms of the mixture to select what we call skeleton terms and the other uses orthogonalization for the same purpose. Both algorithms require O ( r 2 N + p ( d ) r N ) operations, where N is the initial number of terms in the multivariate mixture, r is the number of selected linearly independent terms, and p ( d ) is the cost of computing the inner product between two terms of a mixture in d variables. For general Gaussian mixtures p ( d ) ∼ d 3 since we need to diagonalize a d × d matrix, whereas for separated representations p ( d ) ∼ d (there is no need for diagonalization). Due to conditioning issues, the resulting accuracy is limited to about one half of the available significant digits for both algorithms. We also describe an alternative algorithm that is capable of achieving higher accuracy but is only applicable in low dimensions or to multivariate mixtures in separated form. We describe a number of initial applications of these algorithms to solve partial differential and integral equations and to address several problems in data science. For data science applications in high dimensions, we consider the kernel density estimation (KDE) approach for constructing a probability density function (PDF) of a cloud of points, a far-field kernel summation method and the construction of equivalent sources for non-oscillatory kernels (used in both, computational physics and data science) and, finally, show how to use the new algorithm to produce seeds for subdividing a cloud of points into groups.

We introduce a new functional representation of probability density functions (PDFs) of non-negative random variables via a product of a monomial factor and linear combinations of decaying exponentials with complex exponents. This approximate representation of PDFs is obtained for any finite, user-selected accuracy. Using a fast algorithm involving Hankel matrices, we develop a general numerical method for computing the PDF of the sums, products, or quotients of any number of non-negative random variables yielding the result in the same type of functional representation. We present several examples to demonstrate the accuracy of the approach.

We present a new fast algorithm for computing the Boys function using a nonlinear approximation of the integrand via exponentials. The resulting algorithms evaluate the Boys function with real and complex valued arguments and are competitive with previously developed algorithms for the same purpose., Comment: 9 pages, two figures and two tables

By using Poisson's summation formula, we calculate periodic integrals over Gaussian basis functions by partitioning the lattice summations between the real and reciprocal space, where both sums converge exponentially fast with a large exponent. We demonstrate that the summation can be performed efficiently to calculate two-center Gaussian integrals over various kernels including overlap, kinetic, and Coulomb. The summation in real space is performed using an efficient flavor of the McMurchie-Davidson recurrence relation. The expressions for performing summation in the reciprocal space are also derived and implemented. The algorithm for reciprocal space summation allows us to reuse several terms and leads to a significant improvement in efficiency when highly contracted basis functions with large exponents are used. We find that the resulting algorithm is only between a factor of 5 and 15 slower than that for molecular integrals, indicating the very small number of terms needed in both the real and reciprocal space summations. An outline of the algorithm for calculating three-center Coulomb integrals is also provided.

We report the first fully numerical approach for relativistic quantum chemical calculations applicable to molecules. The approach uses an adaptive basis of multiwavelet functions to solve the full four-component Dirac-Coulomb equation to a user-specified accuracy. The accuracy of the code is demonstrated by comparison with ground state energy calculations of atoms performed in GRASP, and the applicability to molecules is shown via ground state calculations of some simple molecules, including water analogs up to H2Po. In the case of molecules, comparison is made with Gaussian basis set calculations in DIRAC.

We construct high-order derivative operators for smooth functions represented via discontinuous multiwavelet bases. The need for such operators arises in order to avoid artifacts when computing functionals involving high-order derivatives of solutions of integral equations. Previously high-order derivatives had to be formed by repeated application of a first-derivative operator that, while uniquely defined, has a spectral norm that grows quadratically with polynomial order and, hence, greatly amplifies numerical noise (truncation error) in the multiwavelet computation. The new constructions proceed via least-squares projection onto smooth bases and provide substantially improved numerical properties as well as permitting direct construction of high-order derivatives. We employ either b-splines or bandlimited exponentials as the intermediate smooth basis, with the former maintaining the concept of approximation order while the latter preserves the pure imaginary spectrum of the first-derivative operator and provides more direct control over the bandlimit and accuracy of computation. We demonstrate the properties of these new operators via several numerical tests as well as application to a problem in nuclear physics. Keywords: Multiwavelets, Multiresolution, Derivatives, Numerical, Discontinuous, Bandlimited exponentials

We present a new adaptive method for electronic structure calculations based on novel fast algorithms for reduction of multivariate mixtures. In our calculations, spatial orbitals are maintained as Gaussian mixtures whose terms are selected in the process of solving equations. Using a fixed basis leads to the so-called basis error since orbitals may not lie entirely within the linear span of the basis. To avoid such an error, multiresolution bases are used in adaptive algorithms so that basis functions are selected from a fixed collection of functions, large enough to approximate solutions within any user-selected accuracy. Our new method achieves adaptivity without using a multiresolution basis. Instead, as part of an iteration to solve nonlinear equations, our algorithm selects the ‘best’ subset of linearly independent terms of a Gaussian mixture from a collection that is much larger than any possible basis since the locations and shapes of the Gaussian terms are not fixed in advance. Approximating an orbital within a given accuracy, our algorithm yields significantly fewer terms than methods using multiresolution bases. We demonstrate our approach by solving the Hartree–Fock equations for two diatomic molecules, HeH + and LiH, matching the accuracy previously obtained using multiwavelet bases.

We introduce a new method for functional representation of oscillatory integrals within any user-supplied accuracy. Our approach is based on robust methods for nonlinear approximation of functions via exponentials. The complexity of evaluation of the resulting representations of the oscillatory integrals does not depend or depends only mildly on the size of the parameter responsible for the oscillatory behavior.

We present for the first time real-space, arbitrarily-accurate representations of the operators required for up to second-order Douglas-Kroll-Hess (DKH), a model for constructing quasi-relativistic electronic Hamiltonians. The approach can be extended to other operator-based quasi-relativistic models. The representations are in the form of sums of Gaussian functions with positive coefficients and thus enable efficient and numerically-accurate formulations using conventional Gaussian basis sets or other bases such as multiwavelets. The operators are demonstrated with application to hydrogen-like systems using the relativistic-kinematic and first-order DKH Hamiltonians.

The standard particle-in-cell algorithm suffers from grid heating. There exists a gridless alternative which bypasses the deposition step and calculates each Fourier mode of the charge density directly from the particle positions. We show that a gridless method can be computed efficiently through the use of an Unequally Spaced Fast Fourier Transform (USFFT) algorithm. After a spectral field solve, the forces on the particles are calculated via the inverse USFFT (a rapid solution of an approximate linear system). We provide one and two dimensional implementations of this algorithm with an asymptotic runtime of $O(N_p + N_m^D \log N_m^D)$ for each iteration, identical to the standard PIC algorithm (where $N_p$ is the number of particles and $N_m$ is the number of Fourier modes, and $D$ is the spatial dimensionality of the problem) We demonstrate superior energy conservation and reduced noise, as well as convergence of the energy conservation at small time steps., 17 pages, 12 figures

The standard particle-in-cell algorithm suffers from finite grid effects which break energy conservation, cause numerical dispersion, and create numerical instabilities. We present a gridless alternative, bypassing the deposition step and calculating each Fourier mode of the charge density directly from the particle positions. This can be done efficiently through the use of an Unequally Spaced Fast Fourier Transform (USFFT) algorithm1 2. After a spectral field solve, the forces on the particles are calculated via the inverse USFFT (a rapid solution of an approximate linear system). We provide a 1D implementation of this algorithm with an asymptotic runtime of O(Np + Nm log Nm) for each iteration, identical to the standard PIC algorithm (where Np is the number of particles and Nm is the number of Fourier modes). We demonstrate superior energy conservation and reduced noise, as well as convergence at small time steps.

A fully numerical method for the time-dependent Hartree-Fock and density functional theory (TD-HF/DFT) with the Tamm-Dancoff (TD) approximation is presented in a multiresolution analysis (MRA) approach. From a reformulation with effective use of the density matrix operator, we obtain a general form of the HF/DFT linear response equation in the first quantization formalism. It can be readily rewritten as an integral equation with the bound-state Helmholtz (BSH) kernel for the Green's function. The MRA implementation of the resultant equation permits excited state calculations without virtual orbitals. The integral equation is efficiently and adaptively solved using a numerical multiresolution solver with multiwavelet bases. Our implementation of the TD-HF/DFT methods is applied for calculating the excitation energies of H2, Be, N2, H2O, and C2H4 molecules. The numerical errors of the calculated excitation energies converge in proportion to the residuals of the equation in the molecular orbitals and response functions. The energies of the excited states at a variety of length scales ranging from short-range valence excitations to long-range Rydberg-type ones are consistently accurate. It is shown that the multiresolution calculations yield the correct exponential asymptotic tails for the response functions, whereas those computed with Gaussian basis functions are too diffuse or decay too rapidly. We introduce a simple asymptotic correction to the local spin-density approximation (LSDA) so that in the TDDFT calculations, the excited states are correctly bound.

We introduce a new approximate multiresolution analysis (MRA) using a single Gaussian as the scaling function, which we call Gaussian MRA (GMRA). As an initial application, we employ this new tool to accurately and efficiently compute the probability density function (PDF) of the product of independent random variables. In contrast with Monte-Carlo (MC) type methods (the only other universal approach known to address this problem), our method not only achieves accuracies beyond the reach of MC but also produces a PDF expressed as a Gaussian mixture, thus allowing for further efficient computations. We also show that an exact MRA corresponding to our GMRA can be constructed for a matching user-selected accuracy.

We introduce a new, quadratically convergent algorithm for finding maximum absolute value entries of tensors represented in the canonical format. The computational complexity of the algorithm is linear in the dimension of the tensor. We show how to use this algorithm to find global maxima of non-convex multivariate functions in separated form. We demonstrate the performance of the new algorithms on several examples., Comment: 13 pages, 4 figures

This paper introduces a randomized variation of the alternating least squares (ALS) algorithm for rank reduction of canonical tensor formats. The aim is to address the potential numerical ill-conditioning of least squares matrices at each ALS iteration. The proposed algorithm, dubbed randomized ALS, mitigates large condition numbers via projections onto random tensors, a technique inspired by well-established randomized projection methods for solving overdetermined least squares problems in a matrix setting. A probabilistic bound on the condition numbers of the randomized ALS matrices is provided, demonstrating reductions relative to their standard counterparts. Additionally, results are provided that guarantee comparable accuracy of the randomized ALS solution at each iteration. The performance of the randomized algorithm is studied with three examples, including manufactured tensors and an elliptic PDE with random inputs. In particular, for the latter, tests illustrate not only improvements in condition numbers, but also improved accuracy of the iterative solver for the PDE solution represented in a canonical tensor format.

SUMMARY We present a survey sinking based framework for residual migration velocity analysis. By using a migration algorithm in the temporal frequency domain, our migration algorithm is efficient and allows us to pick velocity perturbations directly. Picking velocity perturbation directly removes the issues associated with estimating velocity perturbations from indirect picking parameters, such as time lags. Our methodology is based on a wave-equation migration method that combines high accuracy with computational efficiency. By relying on survey sinking migration, we avoid the problem of estimating the source wavelet, which is a rarely mentioned problem associated with RTM and shot-gather migration. Our survey sinking framework allows us to robustly estimate the velocity sequentially in relatively thin layers.

MADNESS (multiresolution adaptive numerical environment for scientific simulation) is a high-level software environment for solving integral and differential equations in many dimensions that uses adaptive and fast harmonic analysis methods with guaranteed precision that are based on multiresolution analysis and separated representations. Underpinning the numerical capabilities is a powerful petascale parallel programming environment that aims to increase both programmer productivity and code scalability. This paper describes the features and capabilities of MADNESS and briefly discusses some current applications in chemistry and several areas of physics.