The Pseudospectral Method

Like conventional ab initio electronic structure codes, Jaguar solves the Schrödinger equation iteratively, using self-consistent field methods to calculate the lowest-energy wave function within the space spanned by the selected basis set. For calculations on large molecules, both conventional and pseudospectral techniques must recalculate key integral terms for each SCF iteration, since storage costs for these terms are prohibitive.

Most of the fundamental integrals calculated in the pseudospectral method [1‑9] are computed in physical space, on a grid, rather than in the spectral space defined by the basis functions. The pseudospectral method takes the density matrix from the wave function at the beginning of each SCF iteration and the values of the integrals on the grid points and manipulates them to produce the necessary operators on the grid, then assembles the Fock matrix by transforming these components back into spectral space, where the Fock matrix is used in the usual way to generate the wave function for the next iteration.

For medium and large molecules, the additional overhead of the pseudospectral method in computing the transformation between physical and spectral space is vastly outweighed by the advantages of evaluating the integrals in physical space. The matrix needed for the transformation from physical to spectral space [7] can be assembled before the SCF iterations by calculating the least-squares operator Q, which is given by the equation

(1)

where S is the analytic overlap matrix between the fitting functions and the basis set, R is the matrix of fitting functions evaluated at the grid points, and w is a diagonal matrix of grid weights. The fitting functions used to construct the matrix R include both basis functions and dealiasing functions, which are chosen in order to span the function space represented by the grid more completely than the basis functions alone. The operator Q can be calculated for the relevant basis functions using several different sets of grid points, where each set of points defines a grid type, ranging from coarse to ultrafine.

In practice, not all possible elements are calculated for each basis function i and each grid point g, because most basis functions drop off sharply enough that they have no significant value on some or most grid points. These functions are classified as short-range functions and are grouped together by atom, while the remaining functions are classified as long-range functions, which are all considered to be in one single group [13].

Since Q does not depend on the wave function itself, it can be fully computed before the SCF procedure. However, since the Q for each grid type contains N_basis x N_grid elements, where N_basis is the number of basis functions and N_grid the number of grid points (which is generally larger than N_basis), we sometimes reduce memory demands by only computing and storing the N_basis x N_fit matrix in the program rwr, for cases where the Q for that grid type is only needed for one SCF iteration. We then assemble the full Q during the SCF iteration for which it is needed.

After the program rwr has generated the Q or matrix, the program scf takes the initial orbitals and iteratively modifies them with the pseudospectral method until convergence. This process involves calculating the values of the necessary integrals on the grid points, and actually assembling the Fock matrix from the computed information. The three-center, one-electron pseudospectral integrals on the grid points are defined by

(2)

where and are basis functions and the index g represents a grid point. These integrals are calculated for all combinations of basis functions and grid points not eliminated by cutoffs, and the Fock matrix is assembled from its Coulomb and exchange matrix components and , which are calculated in physical space and transformed back into spectral space by the following equations:

(3)

(4)

where D is the usual spectral space density matrix, is the value of the function j at grid point g, and is given by Eq. (2). The grid points used for each SCF iteration are determined by the grid type (coarse, medium, fine, or ultrafine) chosen for that iteration. The number of arithmetic operations involved in the assembly of the matrices J and K in Eq. (3) and Eq. (4) scales formally as N³, as opposed to the N⁴ scaling for the matrix assembly in the conventional spectral space algorithm.

Jaguar actually uses the pseudospectral method described above for the majority of the computationally intensive two-electron integral terms, but calculates the one-electron and some of the largest and most efficiently computed two-electron terms analytically [13]. For the Coulomb matrix elements, we calculate the analytic terms

for cases in which i, j, k, and l meet certain cutoff criteria and the two-electron integral (ij|kl) is of the form (aa|aa), (aa|ab), (aa|bb), (ab|ab), or (aa|bc), where a, b, and c indicate the atom upon which the function is centered. Similar correction terms are computed for the exchange operator, as detailed in ref. [13]. The corresponding pseudospectral terms, as defined by Eq. (3) and Eq. (4) for the appropriate choices of i, j, k, and l, must be subtracted from the pseudospectral J and K elements as well. This combined pseudospectral/analytic approach allows Jaguar to take advantage of the strengths of both methods, since it can largely maintain the pseudospectral method speedups for a particular grid, and can also use a coarser grid than a purely numerical calculation would allow.