System-Specific Separable Basis Based on Tucker Decomposition: Application to Density Functional Calculations.

For fast density functional calculations, a suitable basis that can accurately represent the orbitals within a reasonable number of dimensions is essential. Here, we propose a new type of basis constructed from Tucker decomposition of a finite-difference (FD) Hamiltonian matrix, which is intended to reflect the system information implied in the Hamiltonian matrix and satisfies orthonormality and separability conditions. By introducing the system-specific separable basis, the computation time for FD density functional calculations for seven two- and three-dimensional periodic systems was reduced by a factor of 2-71 times, while the errors in both the atomization energy per atom and the band gap were limited to less than 0.1 eV. The accuracy and speed of the density functional calculations with the proposed basis can be systematically controlled by adjusting the rank size of Tucker decomposition.

Introduction

Numerical methods for replacing partial differential equations with finite-dimensional algebraic equations are essential for the electronic structure calculations of molecular or solid systems.[1,2] The atom-centered and plane-wave basis methods are frequently used for nonperiodic and periodic systems, respectively.[3,4] Real-space methods are potentially competitive with the aforementioned methods because of their flexibility and computational simplicity.[5] However, they have not yet been widely adopted for electronic structure calculations. Discretization of the simulation domain results in a large dimension for the Hamiltonian and orbitals. To mitigate the increase in memory usage and computational costs due to the large dimension, tensor decomposition techniques can be applied to real-space methods.[6−9]

Tensor decomposition techniques are not limited to real-space methods.[10] They have been actively investigated to accelerate various quantum chemistry methods that require a large amount of computational and memory resources (e.g., a perturbation method,[10,11] coupled cluster theory,[12−14] and full basis representation methods[15]). For density functional calculations, tensor decomposition techniques using real-space methods have been studied because orbital values on a rectangular grid can be represented as an order-3 tensor. Solala et al. applied tensor decomposition to density functional calculations to minimize memory load.[6] Their results show that the Tucker decomposition method can successfully compress the orbitals represented on a cubic grid from the results of bubbles and the cube numerical framework, which is a variation of real-space methods. Tensor decomposition can be used to compress orbitals on a three-dimensional (3D) grid and build a basis for a self-consistent field (SCF) procedure. Gavini et al. proposed a Tucker tensor basis derived from a separable approximation of the Hamiltonian. It effectively reduces the dimensions of a Kohn–Sham (KS) Hamiltonian matrix originally represented on an equidistant finite-element grid.[8]

The separability of the basis is an important property to reduce the computational costs of many operations. A typical separable basis can be obtained from the simple product of three 1D functions (e.g., a Gaussian function). To impose system information on such a basis, a 1D function that reflects the system information on a general polyatomic structure needs to be developed. By contrast, a Hamiltonian matrix that implicitly includes all system information can be easily constructed on a real-space grid. Using the Hamiltonian on a real-space grid is an attainable solution for imposing system information on a separable basis.

Herein, we propose a system-specific separable basis derived from a finite-difference (FD) KS Hamiltonian matrix and investigate its performance for 2D and 3D periodic structures. The resulting basis is constructed by reflecting information on the Hamiltonian of the system and is also separable along the axes of the spatial coordinates. These two features are common to other types of Tucker tensor basis.[8] A key contribution of this work is that the basis is built directly from a finite-difference method, and its nonzero patterns are used in the projection process instead of introducing a separable Hamiltonian. In addition, the convergence of our basis is systematically controlled by increasing the rank size of Tucker decomposition. In the following, we briefly explain the mathematical background of our method, followed by the implementation details. We then discuss the performance of the proposed basis on density functional calculations for 2D and 3D periodic systems and demonstrate its advantages for reducing the computation time of density functional calculations.

Method

Tucker Representation and Higher-Order Singular-Value Decomposition

Here, we briefly introduce a Tucker representation and a higher-order singular-value decomposition (HOSVD) method for completeness. A more detailed explanation can be found in previous papers.[16,17]

The Tucker representation is used to represent an order-d tensor, , as a contraction of a small order-d core tensor, , and d unitary factor matrices, (n ∈ {1, 2, ..., d}), where N and r are the dimensions of the nth axis for the original and core tensors, respectively.[16−19] Then, the Tucker decomposition can be written asThe convergence of eq is mathematically guaranteed as r approaches N.[10] However, the Tucker representation is frequently used to find a compact representation of a given tensor, which means r < N. This compact representation can reduce the computational complexity and memory consumption of tensor operations, thereby minimizing accuracy loss.[17]

An HOSVD method is the most common choice to find a set of U( and the corresponding .[10,16] This is one multilinear extension of the matrix singular-value decomposition (SVD). In the HOSVD method, U( is obtained as the left singular vector of a factor-n flattened tensor . Because all U( from HOSVD are unitary, in eq can be evaluated from the contraction of the original tensor with the obtained U(, as follows:where (·)* indicates a complex conjugate.

To obtain a compact Tucker representation, the singular vectors of with large singular values are denoted as . The compact core tensor is computed using consecutive tensor contractions, as shown in eq . Although the memory usage for and U( is much smaller than that for , the key patterns of can be recovered by contraction with U(, as in a typical compact SVD.

In the field of quantum chemistry, the Tucker representation has been used to accelerate the calculation of two-electron integrals of atom-centered basis functions or tensor contractions for higher-order methods.[10,12,20−22] Here, we build a separable basis using the Tucker decomposition of the KS Hamiltonian in the FD representation.

Tucker Decomposition of the Finite-Difference Hamiltonian

For an N × N × N Cartesian grid, the FD Hamiltonian matrix, H, is constructed on . H can be reshaped into an order-6 tensor , the elements of which are given bywhere p and q denote the indices of grid points whose x, y, and z indices are (i, j, k) and (i′, j′, k′), respectively.

To ensure that the factor matrices of the Hamiltonian matrix span low-lying orbitals well, we introduce a constant fictitious potential, which implies that . This constant potential shifts the eigenspectrum of the original Hamiltonian downward without changing the eigenvectors, so that the factor matrices from the decomposition of a Hamiltonian matrix with fictitious potentials are more likely to span low-lying orbitals of the Hamiltonian, which are physically meaningful. For simplicity, we do not denote the fictitious potential in this section. A more detailed explanation for the fictitious potential is provided in the Appendix.

Using the HOSVD method, can be decomposed into the core tensor, , and the corresponding factor matrices, U:where α, β, γ, α′, β′, and γ′ are the indices of .

Owing to the Hermitian property of H, if r1 = r4, r2 = r5, and r3 = r6, the flattened Hamiltonian matrices and factor matrices satisfy the following relations:Hereafter, for convenience, we use U, U, U, H(, H(, and H( instead of U(1), U(2), U(3), H(1), H(2), and H(3), respectively. Similarly, the rank sizes of U, U, and U are denoted as r, r, and r, respectively.

However, we define the square matrix form of aswhere μ and ν are the indices of . From eqs and 5, H̃ can be rewritten aswhere ⊗ and (·) denote the Kronecker product and the conjugate transpose, respectively. Here, H̃ is the projection of H on separable basis vectors U (≔ U ⊗U ⊗U). As discussed in the previous section, the convergence of U, U, and U to make both sides of eq equal is mathematically guaranteed as r, r, and r reach N, N, and N, respectively. Therefore, it is guaranteed that H̃ becomes equal to H when its dimension r × r × r becomes N × N × N.

Here, U is a set of separable basis vectors that can reduce the dimensions of the Hamiltonian from N × N × N to r × r × r. In addition, U satisfies the orthonormality condition because U, U, and U are orthonormal matrices. If U spans physically meaningful eigenstates of the original Hamiltonian well, we only need to diagonalize H̃, which has a smaller dimension than that of H. In addition, U denotes a set of numerical basis vectors that are never explicitly constructed. Owing to its separability, operations with U can be replaced by operations with three small matrices, U, U, and U. Therefore, the memory requirement for U is not N × N × N × r × r × r but rather N × r + N × r + N × r.

To evaluate H̃, instead of directly evaluating the right-hand side of eq , we project three terms of the Hamiltonian matrix separately: kinetic energy, local potential, and nonlocal potential terms. Owing to the properties of U and the nonzero patterns of the three terms, the evaluation of H̃ can be efficiently performed. A further explanation of the projection of the Hamiltonian matrix is described in the Appendix.

Unlike typical basis functions that use a predetermined formula, our U reflects the system information (e.g., relative positions of atoms, phase factors, and cell size) because it is constructed from the decomposition of the Hamiltonian matrix that includes system information. Hereafter, we name U the system-specific separable basis vector and investigate its applicability to density functional calculations.

Before we discuss the performance of the system-specific separable basis in a density functional calculation, we plot its overall process in Figure . The right and left panels of Figure represent the conventional SCF procedure and the additional process, respectively. In the system-specific basis calculation, the basis vectors are constructed using the eigendecomposition of H(H(, H(H(, and H(H(. The eigenvectors of H(H(, H(H(, and H(H( are identical with the left singular vectors of H(, H(, and H(, respectively. To avoid the SVD of a large sparse matrix, we perform eigendecomposition instead of SVD.

Figure 1

Schematic illustration of a conventional self-consistent field procedure (right panel) and the additional process introduced by the system-specific separable basis (left panel).

After basis construction, the original Hamiltonian is projected to the obtained basis space. The eigenvalues and eigenvectors of the projected Hamiltonian are then computed using a typical matrix diagonalization method. The orbital, ϕ, for computing the density, ρ, is evaluated as , where ϕ and are the pth eigenvectors of H and H̃, respectively. We note that ϕ satisfies the orthonormality condition because both and U are orthonormal. The Hartree and the exchange-correlation (XC) potentials for the obtained ρ are evaluated in the same way as the ordinary FD calculation.

To accelerate the system-specific basis vector calculations, we introduce two approximations. The first approximation is using the fixed basis during the SCF loop. In other words, the basis is constructed in the first step of the SCF loop using the initial Hamiltonian matrix; it is then used in the subsequent SCF steps. Although the construction of the basis set is not computationally heavy, changes in the basis set at every SCF step reduce the speed of SCF convergence. For a fixed basis, only two local potential terms (Hartree and XC potentials) must be updated at each SCF step. Therefore, only the projection of the updated local potential is performed for each SCF step.

The second approximation is discarding the nonlocal pseudopotential in the basis construction. The errors introduced by the two approximations are plotted in Figure S1. The approximations may induce errors up to 100 meV in both the atomization energy and the band gap; however, these deviations disappear when the rank size sufficiently converges, and the approximations lead to a ∼6× increase in speed in all tested cases. Hereafter, all results are obtained using both approximations.

Implementation and Experiments

The construction of system-specific basis vectors and the projection of the Hamiltonian are performed using the Tucy package, which is written in C++ and has a Python interface. For density functional calculations, our Python package, called the grid-based open-source Python engine for large-scale simulation (GOSPEL), was used. GOSPEL supports FD calculations and system-specific separable basis calculations using the Tucy. In GOSPEL, the Hartree potential was evaluated using the interpolation scaling method, as in our previous studies.[23,24] An XC potential is evaluated from the experimental version of libXC.[25]

To assess the convergence and performance enhancement, both the reference FD and system-specific separable basis calculations were performed using the same systems. All computational options were used equally in both cases, and all calculations were performed using a single thread of an Intel Xeon Gold 6234 CPU. The PBE[26] functional was used for the XC functional, and optimized norm-conserving Vanderbilt[27] pseudopotentials were used. All 2D and 3D periodic structures were calculated using (4 × 4 × 1) and (4 × 4 × 4) k-point meshes, respectively. For the kinetic energy matrix, a 7-point FD matrix is used. SCF procedures end when the sum of the occupied band energies converges to less than 10–6 Hartree.

For iterative diagonalization for both typical FD and the projected Hamiltonian matrices, we use LOBPCG functions implemented in scipy,[28] a highly mature and optimized Python package for scientific computing. A compressed sparse row format is used for the FD Hamiltonian instead of a dense matrix format to compute the matrix-vector multiplications. Tucy and GOSPEL are freely available in their online git repositories (https://gitlab.com/jhwoo15/gospel and https://gitlab.com/jhwoo15/tucy, respectively).

Because the cell parameters do not exactly match the multiples of a given grid spacing, the actual grid spacing is set to have the closest value of the given spacing within a small difference (up to 0.1 bohr). Here, we denote the given grid spacing instead of the actual grid spacing in the paper for better readability. The actual grid spacing corresponding to each structure is listed in Table S1.

The performance of the system-specific separable basis is assessed for seven structures: three cubic diamond structures (C, SiC, and Si), two ABO3 perovskites (BaTiO3 and SrTiO3), and two hexagonal 2D materials (hBN and hBCN). The atomic coordinates and cell parameters of the seven systems are presented in the third section of the Supporting Information. We used the atomization energy per atom for a fair comparison between systems with different numbers of atoms. Hereafter, we refer to the atomization energy per atom as the atomization energy.

Results and Discussion

To confirm the dependency of the system-specific separable basis on the system information, we plot the results of the 1D Fourier transform of the first five vectors of the factor matrices constructed from the initial Hamiltonian at the Γ- (blue bars) and X- (red bars) points (see Figure ). In the lowest panel, only Γ-point data is visible because only one k-point was sampled along the nonperiodic axis of the 2D hexagonal sheet structures.

Figure 2

Results of 1D Fourier transform of the first five vectors of factor matrices from the initial Hamiltonian of the seven different systems: C, SiC, Si, SrTiO3, BaTiO3, hBN, and hBCN. The blue and red dashed lines indicate basis vectors from the Γ- and X-points, respectively. U and U denote the factor matrices of the x- and z-axis, respectively, where the z-axis is a nonperiodic axis of hBN and hBCN.

First, the blue bars are symmetrically distributed in all cases because the basis vectors at the Γ-point are always real regardless of the structure. However, at the X-point, the basis vectors and Hamiltonian matrices are no longer symmetric because of the phase factor; therefore, the basis vectors are no longer symmetric in Fourier space. One interesting point related to the k-space is that the basis vectors at the X-point are not just shifts in the basis vector at the Γ-point. This indicates that the basis vectors at different points in the k-space are not the products of Γ-point basis vectors with the phase factor, and the basis vectors at each k-point are constructed in a way that reflects the overall Hamiltonian matrix.

Figure also shows the structural dependency of the basis vectors. The ABO3 structures have a common pattern. C and Si structures also share a similar pattern. However, the SiC structure has a different shape than those of other diamond structures. SiC is composed of two different elements; therefore, the nature of the covalent bonds in SiC is largely different than those of C and Si. Likewise, the basis vectors of two hexagonal sheet structures show different patterns along the x-axis, whereas the basis vectors along the z-axis show a similar trend. This implies that hBN and hBCN show different characteristics along the periodic axes but not along the nonperiodic axis. Although it is difficult to elucidate which system information changes a specific pattern in the basis vectors that we obtain, we can observe the structural and phase dependencies of the system-specific basis vectors.

We investigated whether the obtained basis vectors can properly span an orbital from the reference FD calculations. We projected the reference orbitals from the FD calculation of the SrTiO3 on U and calculated their residuals. Figure plots the sizes of the projection residual on U constructed with different rank sizes (r, r, and r). The large residual size means that the basis space does not sufficiently span the reference orbital. The tested reference orbitals were obtained from ordinary FD calculations of SrTiO3. For good readability, we present the residual size of the first 100 orbitals and occupied orbitals in Figure a and 3b, respectively.

Figure 3

Size of the projection residual of (a) the first 100 orbitals and (b) the occupied orbitals of SrTiO3 at the Γ-point. The residual is computed by the projection of the reference orbitals, ϕ, on the separable basis vectors, U, obtained from Tucker decomposition with different rank sizes.

For small rank sizes, the basis vectors do not sufficiently cover the reference orbitals, but they span the orbitals better as the rank size increases (see Figure a). In addition, it was observed that the residual sizes for the low-lying orbitals do not always decrease first as the rank size increases, but those for high virtual orbitals slowly converge to zero. This indicates that the basis space spans the orbitals sufficiently, especially for low-lying orbitals.

To investigate the effects of system-specific separable basis vectors in the SCF procedure, we performed density functional calculations for the tested structures. The computational details and results are summarized in Tables S3–S9.

Figure shows the absolute errors in atomization energies, |ΔEa| (black line), and band gaps, |ΔEg| (blue line), as a function of the number of basis vectors, R = r × r × r. For the cubic diamond structures and ABO3 structures, we sampled the same rank sizes for all three axes, whereas the rank size of the hexagonal sheet structures is proportional to the cell size along each axis. The same cell parameters and rank sizes were used for the single-atom calculations needed to calculate the atomization energies.

Figure 4

Convergence behavior of system-specific separable basis vectors as a function of the number of basis vectors, R, on a log–log scale. The blue and black lines indicate the absolute errors in the band gap and atomization energy per atom (|ΔEg| and |ΔEa|), respectively. Both the system-specific separable basis and reference finite-difference calculations were performed with h = 0.2 bohr.

The atomization energy and the band gap do not converge monotonically with respect to R, whereas the monotonic convergence of the total energy is guaranteed by the variational principle, as shown in the 10th column in Tables S3–S9. In the test range of R, the systems show convergence within 0.1–10 meV for both |ΔEa| and |ΔEg|. Elucidating the dependence of the error convergence on the structures is difficult with a few test cases. Nonetheless, we note that the system-specific separable basis converges well, even in systems with transition metals or a nonperiodic axis.

To investigate the speed of SCF calculations with a system-specific separable basis, we plot the elapsed times for the overall calculations and three major bottlenecks for Si calculations as a function of R in Figure . To check the dependence of the computation time on the grid spacing h, we also plotted the results with different h values. The blue, orange, and green lines indicate the results with h values of 0.3, 0.25, and 0.2 bohr, respectively. The dashed lines represent the elapsed time for the reference FD calculations. The elapsed time of the total SCF procedure with the system-specific separable basis increases as R increases but does not show a dramatic change with respect to h, as shown in Figure a. By contrast, the total elapsed time of the reference FD results increased rapidly as h decreased.

Figure 5

Elapsed times of density functional calculations according to the number of basis vectors on a log–log scale. Elapsed time for (a) total SCF calculation, (b) diagonalization, (c) projection of the Hamiltonian matrix, and (d) construction of the system-specific separable basis. The dashed lines in (a) and (b) represent the elapsed time for the reference finite-difference calculations for the same system (Si) and computational conditions. Blue, orange, and green lines represent the results for h of 0.3, 0.25, and 0.2 bohr, respectively.

The increase in the elapsed time of the reference calculations originates from diagonalization, which is the primary bottleneck. As shown in Figure b, most of the elapsed time for the reference calculations is spent in diagonalization, and its cost is strongly dependent on the h values. For the case of a system-specific separable basis, the elapsed time of the diagonalization is independent of the choice of h and is much smaller than that of the reference cases. This is because the dimensions of the projected Hamiltonian matrix are determined not by the number of grid points, N = N × N × N, but by R which is much smaller than N.

The system-specific separable basis additionally induces the basis construction and projection processes. Figures c and d show the elapsed time for projection and basis construction, respectively. The computational time for basis construction relies on h values because we obtain U, U, and U from the direct diagonalization of small dense matrices, H(H(,H(H(, and H(H(. Despite the strong h dependence, the basis construction time occupies only a small part of the overall time. Contrary to the basis construction time, the cost of the projection depends on both h and R. The detailed computational complexities of the projections and basis construction are explained in the Appendix. The projection time increases as R increases and h decreases. However, large differences in the elapsed time of the projection as a function of h are not shown, except for a few small R cases. Hence, the total computational time for a system-specific basis calculation does not increase significantly as h decreases.

Although system-specific separable basis calculations require additional processes, they show excellent performance in diagonalization; thus, the overall computational costs are reduced in most cases. Here, we discuss only the results of Si, but we observed the same trend for other systems (see the sixth–ninth columns of Tables S3–S9).

Figure summarizes the overall performance enhancement by the new basis with respect to the reference FD calculations as a function of R/N, when h = 0.2 bohr. Figure S3 shows the performance enhancement results with other h values. For all systems, smaller R/N values resulted in larger performance enhancement. This is because a smaller R/N implies a larger reduction in the dimension of the Hamiltonian matrix, resulting in faster diagonalization. The elapsed time for each component for all calculations is included in Tables S3–S9. The intersections of the horizontal and vertical lines represent the smallest R/N case for each system, where both |ΔEa| and |ΔEg| were less than 100 meV. The intersections of the ABO3 and hexagonal sheet structures were ∼8% and ∼2%, respectively. The diamond structures should have an R/N of ∼5% to achieve a tolerance of 100 meV, except for Si. The N value for the Si system is greater than those of the other diamond structures because the cell volume of Si is larger than that of the others. In addition, the R for the error convergences was slightly smaller than that of the others. Therefore, the Si system showed significant performance enhancement.

Figure 6

Performance enhancement by a system-specific separable basis as a function of the ratio of dimensions of the finite-difference and the projected Hamiltonian matrices, R/N. The intersections of the horizontal and vertical lines represent the smallest R/N cases, where both the atomization energy and the band gap have errors of less than 100 meV. For all calculations, with or without basis vectors, a grid spacing of 0.2 bohr was used.

To be more practical, a system-specific separable basis must achieve performance improvements with sufficiently high accuracy. The system-specific separable basis balances speed and accuracy by tuning R. A large R achieves high accuracy but simultaneously reduces the calculation speed, as shown in Figures and 6. Table summarizes the performance enhancements of the tested systems for 100, 50, and 25 meV tolerances for both types of errors. As shown in Figure , the use of large R/N to achieve high accuracy reduces the gains in computation speed. However, significant acceleration (2–14×) was achieved, even within a small tolerance value of 25 meV.

Table 1

Performance Enhancements of the Tested Systems at Various Error Tolerancesa

	diamond			ABO3		hexagonal sheet
tolerance (meV)	C	SiC	Si	SrTiO3	BaTiO3	hBN	hBCN
100	5.88 (5.59%)	5.86 (4.56%)	71.43 (0.95%)	3.65 (6.66%)	2.07 (8.95%)	13.31 (2.34%)	24.17 (1.86%)
50	3.14 (6.98%)	2.51 (6.63%)	71.43 (0.95%)	1.65 (9.70%)	2.07 (8.95%)	4.27 (4.53%)	7.70 (3.06%)
25	2.39 (8.59%)	1.71 (7.87%)	14.07 (2.40%)	1.65 (9.70%)	2.07 (8.95%)	4.27 (4.53%)	5.36 (3.79%)

The percentages in parentheses indicate the ratio of dimensions of the finite-difference and the projected Hamiltonian matrices, R/N.

Conclusion

Here, we proposed a system-specific separable basis derived from Tucker decomposition of a finite-difference Hamiltonian and investigate its performance in density functional calculations. We show that the new basis can successfully span low-lying orbitals and that their coverage is systematically improved by increasing the rank size of Tucker decomposition. The proposed basis dramatically reduces the dimensions of the Hamiltonian matrix and hence accelerates the diagonalization of the Hamiltonian matrix. We confirmed the properties of the basis vectors and measured the performance enhancements using seven selected systems. The system-specific separable basis achieved a 2–71× increase in computation speed with 100 meV tolerance for the band gap and the atomization energy. Higher accuracy can be achieved for all tested systems with a larger rank size but a lower gain in computation speed (e.g., 2–14× increase with a 25 meV tolerance). Here, we validated the performance of the system-specific separable basis only for density functional calculations. However, we expect it to be useful for higher-order quantum chemical methods because our basis benefits from the advantages of both the real-space method (e.g., fast numerical integrations and derivatives) and the typical basis function expansion method (e.g., low dimension of basis space and separability).