Hybrid functionals currently used in the generalized KS framework contain a fraction of non-local, exact HFX term. In the PBE0 hybrid functional (also known as PBEh or PBE1PBE) (Adamo and Barone, 1999; Ernzerhof and Scuseria, 1999), the exchange-correlation energy is written as E xc PBE0 = 1 4 E x HF + 3 4 E x PBE + E c PBE (1) where 25% HFX E x HF is mixed with 75% Perdew-Burke–Ernzerhof (PBE) exchange E x PBE , and the electronic correlation is still represented by the part of the PBE correlation E c PBE . The inclusion of the HFX in PBE0 reduces the self-interaction error of the density functional, resulting in a substantial improvement over the parent PBE. However, the full-range (FR) HFX is computationally very demanding and may be problematic in solids. To address this issue, Heyd, Scuseria, and Ernzerhof (Heyd et al., 2003; Heyd et al., 2006) proposed to replace the long-range part of HFX in the PBE0 by a corresponding PBE counterpart. Then, the resulting expression for the HSE exchange–correlation energy is given by E xc HSE = 1 4 E x SR,HF ω + 3 4 E x SR,PBE ω + E x LR,PBE ω + E c PBE (2) where E x SR,HF and E x LR,PBE ( ω ) is the short-range (SR) HFX and long-range (LR) PBE exchange energies, and E x SR,PBE ( ω ) is the SR exchange energy. ω is an adjustable screening parameter that defines the range separation, and ω is set to 0.11 Bohr⁻¹ for HSE06 (Heyd et al., 2006). Such a treatment not only improves the computational convenience but also avoids the problematic effects of LR HFX in metals and semiconductors with narrow band gaps.

For periodic systems, the HFX energy per unit cell can be written as E x HF = − 1 2 ∑ σ = α , β ∑ k , q BZ × ∑ i , j occ ∫ Ω ∫ ψ i k σ * r ψ j q σ r v ̂ r , r ′ ψ j q σ * r ′ ψ i k σ r ′ d r ′ d r (3) where ψ i k σ ( r ) denote the i-th occupied (occ) crystalline spin-orbitals with spin σ for k point sampling in the Brillouin zone (BZ), Ω is the unit cell volume. E x HF is used to represent the FR or SR HFX energy, and v ̂ ( r , r ′ ) is either the Coulomb operator v ̂ ( r , r ′ ) = 1 / | r − r ′ | in PBE0 or the screened Coulomb operator v ̂ ( r , r ′ ) = e r f c ( ω | r − r ′ | ) / | r − r ′ | in HSE. Hereafter, we will formulate only the collinear spin-polarized case with σ = {α, β}.

In the linear combination of atomic orbitals (LCAO) method, the spin-orbitals are expanded in terms of a linear combination of Bloch basis functions ψ i k σ r = 1 N ∑ n N e i k R ∑ μ N b c μ σ k ϕ μ r − r μ − R (4) where ϕ _μ (r − r _μ − R) denotes the μ-th NAO centering at r _μ within a lattice translation vector R, c μ σ ( k ) is the expansion coefficient, and N is the number of primitive unit cells under the Born-von Kármán (BvK) periodic boundary conditions. Within NAOs, the HFX matrix for the self-consistent field (SCF) calculations can be written as H x σ , H F μ κ G = − ∑ ν λ ∑ N , H P ν λ σ , H − N μ 0 ν N | κ G λ H (5) and the corresponding HFX energy can be obtained by E x HF = − 1 2 ∑ σ ∑ μ ν κ λ ∑ G , N , H P μ κ σ , G P ν λ σ , H − N μ 0 ν N | κ G λ H (6) where the subscripts of Greek letters {μ, ν, κ, λ} label NAOs, the superscript 0 represents the reference primitive unit cell, while G, N, and H represent extended unit cells in the BvK supercells. The P μ κ σ , G is the spin density matrix element, which can be obtained by an integration of the expanded coefficients in the BZ P μ κ σ , G = ∑ j ∫ BZ c μ , j σ ∗ k c κ , j σ k θ ϵ F − ϵ j σ k e i k G d k (7) where θ represent the step function, ϵ _F is the fermi energy and ϵ j σ ( k ) is the j-th eigenvalue at k. For hybrid functional calculations with NAOs, the main bottleneck is the evaluation of ERIs μ 0 ν N | κ G λ H = ∬ ϕ μ 0 r ϕ ν N r v ̂ r , r ′ ϕ κ G r ′ ϕ λ H r ′ d r d r ′ (8)

Since the NAOs are dependent of atomic positions, in hybrid functional calculations we must additionally calculate the HFX contribution to atomic forces. TheHFX forces acting on the I-th atom can be directly obtained from the negative gradients of the HFX energy with respect to atomic position R _I F I HF = − ∂ E x HF ∂ R I = ∑ σ ∑ μ κ ∂ P μ κ σ , G ∂ R I ∑ ν λ ∑ N , H P ν λ σ , H − N μ 0 ν N | κ G λ H + 1 2 ∑ σ ∑ μ ν κ λ ∑ G , N , H P μ κ σ , G P ν λ σ , H − N ∂ μ 0 ν N | κ G λ H ∂ R I (9) Note that the first term in Eq. 9 can be rewritten as ∑ μ κ ∂ P μ κ σ , G ∂ R I ∑ ν λ ∑ N , H P ν λ σ , H − N μ 0 ν N | κ G λ H = − ∑ μ κ H x σ , H F μ κ G ∂ P μ κ σ , G ∂ R I (10) which is automatically included in the orthogonalization force due to the non-orthonormality of the NAO basis set (Soler et al., 2002; Li et al., 2016). For periodic systems, the orthogonalization force is given by F I orth = − ∑ σ ∑ μ κ H μ κ σ , G ∂ P μ κ σ , G ∂ R I = ∑ σ ∑ μ κ E μ κ σ , G ∂ S μ κ G ∂ R I . (11) where S μ λ G = ⟨ ϕ μ 0 | ϕ λ G ⟩ is the overlap matrix element, and E μ λ σ , G is the energy-density matrix element given by E μ κ σ , G = ∑ j ∫ BZ c μ , j σ * k c κ , j σ k ϵ j σ k e i k G d k (12) where ϵ j σ ( k ) is the eigenstate energy.

Thus, we only need to deal with the second term 1 2 ∑ σ ∑ μ ν κ λ ∑ G , N , H P μ κ σ , G P ν λ σ , H − N ∂ μ 0 ν N | κ G λ H ∂ R I = 1 2 ∑ σ ∑ μ ν κ λ ∑ G , N , H P μ κ σ , G P ν λ σ , H − N × ∂ ϕ μ 0 ∂ R I ϕ ν N | ϕ κ G ϕ λ H + ϕ μ 0 ∂ ϕ ν N ∂ R I | ϕ κ G ϕ λ H + ϕ μ 0 ϕ ν N | ∂ ϕ κ G ∂ R I ϕ λ H + ϕ μ 0 ϕ ν N | ϕ κ G ∂ ϕ λ H ∂ R I (13) in which the ERI derivatives have to be evaluated properly. Since each ERI may have four different centers, a maximum of 12 differentials is required. Therefore, the calculation of HFX forces is formally more troublesome than that of HFX matrix, and a poor implementation will decrease the overall performance.

A normalized NAO for atom I located at R _I is defined as the product of a numerical radial function and a real regular solid harmonic (Soler et al., 2002) ϕ Ilmζ NAO r = ϕ I l ζ NAO r I r I l Y l m θ , φ (14) where r _I = r − R _I , r _I = |r _I |, l and m label the angular and magnetic momentum quantum numbers, respectively. In multiple-ζ bases, ζ labels different basis with the same quantum numbers (l, m) but different radial shapes. For simplicity, we will omit the index ζ later. The numerical radial function involves a normalization factor N (l, α) and is numerically tabulated in a linear radial mesh. [ r I l Y l m ( θ , φ ) ] also includes its individual normalization factor. NAOs are strictly localized in real space, which provides greater convenience for linear-scaling DFT calculations. However, the evaluation of ERIs over NAOs in real space is computationally expensive, which will introduce a big prefactor in linear-scaling HFX calculations (Shang et al., 2010).

In the LCAO framework, GTOs are by far the most commonly used basis functions to represent the molecular orbitals. This preference is mainly due to the analytical properties of GTOs, which allow efficient evaluation of ERIs for (post-)Hartree-Fock calculations. A normalized primitive spherical harmonic GTO is defined as ϕ Ilmα GTO r = N l , α exp − α r I 2 r I l Y l m θ , φ (15) where α is the orbital exponent, and N (l, α) is the normalization factor over the radial coordinates N l , α = 2 2 l + 3 l + 1 ! α l + 3 / 2 2 l + 2 ! π 1 / 2 1 / 2 (16) It is worth noting that most efficient algorithms for the evaluation of ERIs are based on Cartesian primitive GTOs G l α r = N l x , l y , l z , α x − R x l x x − R y l y z − R z l z ⁡ exp − α r − R I 2 (17) where l = l _x + l _y + l _z labels the angular-momentum quantum number, R _I = (R _x , R _y , R _z ) is the orbital center, and N (l _x , l _y , l _z , α) is the normalization factor N l x , l y , l z , α = 2 π 3 / 4 2 l α 2 l + 3 / 4 2 l x − 1 ‼ 2 l y − 1 ‼ 2 l z − 1 ‼ 1 / 2 (18) For a shell of angular momentum l, there will be 2l + 1 spherical GTOs, but (l + 1) (l + 2)/2 Cartesian GTOs. The transformation between normalized spherical and Cartesian GTOs is required with a transformation matrix c (l, m, l _x , l _y , l _z ) given by Schlegel and Frisch (Schlegel and Frisch, 1995).

In order to obtain ERIs efficiently, we can represent the NAOs in terms of a linear combination of spherical primitive GTOs, and then calculate the ERIs analytically by calling available libraries for GTO-based integrals (e.g. LIBINT (Valeev and Fermann, 2014)). This scheme is called NAO2GTO, which in principle is quite similar to the minimal STO-nG basis set used in the quantum chemistry community. In the NAO2GTO scheme, the numerical radial function of each NAO is fitted as a linear combination of different Gaussians, ϕ l NAO r ≃ ϕ l CGTO r = ∑ i = 1 M D i ⁡ exp − α i r 2 (19) where M is the number of Gaussians, α _i and D _i are the contraction coefficient and exponent, respectively, similar to contracted GTOs (CGTOs).

If the NAO2GTO fitting is strictly accurate, the fitted orbitals will be automatically normalized because of the normalization of NAOs. Once there is a non-negligible fitting error, the fitted orbitals will not be normalized. Since normalization factors between NAOs and CGTOs are different, the following approximation will no longer hold ϕ l m NAO r ≠ ϕ l m CGTO = N l , α , D ∑ i = 1 M D i ⁡ exp − α i r 2 r l Y l m θ , φ (20) with a normalization factor N l , α , D = 2 l + 2 ! π 1 / 2 2 2 l + 3 l + 1 ! ∑ i , j = 1 M D i D j α i + α j l + 3 / 2 − 1 / 2 (21) where α = {α _i } and D = {D _i }. In practice, it is difficult to achieve an exact NAO2GTO fitting with a small number of GTOs, such as M = 3–6. Then, the NAO2GTO fitting will inevitably introduce the ERI errors between original NAOs and fitted CGTOs. In some cases, such errors can even invalidate the final results. In our previous work, we employed the fitted CGTOs to approximately evaluate the ERIs for the HFX term while retaining original NAOs for pure DFT parts. For most systems, we have found that a self-consistent convergence problem often arises in hybrid functional calculations even with high-precision fitting.

To eliminate the fitting errors properly, here we replace original NAOs with the fitted and discretized CGTOs, which is done as the following steps:

• Perform a less rigorous NAO2GTO fitting for each NAO to obtain a set of CGTOs according to Eq. 19;

Note that the fitted CGTOs will give larger cutoff radii compared to the original NAOs. Therefore, we need to feed the values inside a new cutoff radius back to the radial function, beyond which all values are equal to 0.

With the auxiliary CGTOs, one shell set of contracted ERIs (ab|cd) can be calculated by a b | c d = ∑ k K ∑ l L ∑ m M ∑ n N C a k C b l C c m C d n a k b l | c m d n (22) with C a k = D a k c a , m , a x , a y , a z (23) where a denotes the shell orbital with angular momentum a = a _x + a _y + a _z , k is the index of K GTOs, and [a _k b _l |c _m d _n ] represents a set of primitive ERIs over primitive Cartesian GTOs, in which (2a + 1) (2b + 1) (2c + 1) (2d + 1) primitive and contacted ERIs are calculated at once.

Since a primitive ERI contain four centers of A, B, C, and D, its first-order derivatives should have the following 12 terms ∂ a b | c d ∂ A i , ∂ a b | c d ∂ B i , ∂ a b | c d ∂ C i , ∂ a b | c d ∂ D i i ∈ x , y , z . (24) but only 9 derivatives are required because of the translational invariance ∂ a b | c d ∂ A i + ∂ a b | c d ∂ B i + ∂ a b | c d ∂ C i + ∂ a b | c d ∂ D i = 0 (25) Analogously, the first-order derivatives of contracted ERIs can also be evaluated from the primitive ones, which are actually a linear combination of higher and lower angular momentum ERIs ∂ ∂ A i a b | c d = 2 α a + 1 i b | c d − a i a − 1 i b | c d (26)

We obtain the primitive ERIs and their derivatives from the LIBINT library, (Valeev and Fermann, 2014) which implements efficient recursive schemes based on the Obara-Saika method (Obara and Saika, 1986) together with the Head-Gordon-Pople (Head-Gordon and Pople, 1988) and Hamilton-Lindh (Hamilton and Schaefer, 1991; Lindh et al., 1991) variations. We also consider the eight-fold permutational symmetry of ERIs with NAOs, which for periodic systems is given by μ 0 ν N | κ G λ H = μ 0 ν N | λ H κ G = ν 0 μ − N | κ G − N λ H − N = ν 0 μ − N | λ H − N κ G − N = κ 0 λ N − G | μ − G ν N − G = κ 0 λ N − G | ν N − G μ − G = λ 0 κ G − H | μ − H ν N − H = λ 0 κ G − H | ν N − H μ − H (27) In this way, we only need to handle about 1/8 of ERIs and their derivatives. Furthermore, the translational invariance also gives ∂ ( μ μ | μ μ ) ∂ R μ = 0 , which means that we can ignore the ERI derivatives if four orbitals have the same center.

In practice calculations, most of the ERIs and their derivatives have no significant contributions to the HFX matrix and forces, which can be omitted by using integral screening techniques. Thus, integral screening is essential for reducing the computational cost, which should be able to provide an easy-to-estimate upper bound for ERIs. We have previously employed several ERI screening techniques to obtain an efficient and linear-scaling HFX calculation (Shang et al., 2011; Qin et al., 2015). For the calculation of HFX forces, however, integral screening based on the upper bound of ERI derivatives is not a good choice. The reason is that, according to the Schwarz inequality (Häser and Ahlrichs, 1989), the upper bound of ( ∂ μ ∂ R I ν | λ σ ) requires a relatively expensive calculation of ( ∂ μ ∂ R I ν | ∂ μ ∂ R I ν ) (Horn et al., 1991), which does not appear in the ERI derivatives. Alternatively, we can use the same screening techniques based on the upper bound of ERIs as done in the construction of the HFX matrix. That is, if an ERI is skipped during the calculation of the HFX energy, its derivatives should also be neglected for the calculation of HFX forces. Here, we describe all the screening techniques we have employed to compute the HFX forces. For simplicity, we omit the superscripts of 0, G, N, and H later, and the indices {μ, ν, κ, λ} label shell orbitals in the following.

With the rapid development of computer clusters and supercomputers, high-performance computing (HPC) is now essential for program design. In the parallel implementation of HFX forces, the key is to distribute the calculation of ERI derivatives across different central processing unit (CPU) cores. Of course, a straightforward way is to evenly distribute all shell quartets on each processor. However, the total amount of shell quartets is unknown until the integral screening is finished. Furthermore, the computational cost for different types of shell quartets may be very different. For instance, (dd|dd) with higher order angular momentum may be hundreds of times more expensive than (ss|ss), which also may cause serious load imbalance. On the other hand, the distribution of the density matrix may also introduce additional communication when considering the full ERI symmetry. Therefore, load balancing and minimizing communication overhead are essential considerations in the parallel design of HFX force calculation.

For massively parallel computing, a better choice is to employ the master-worker dynamic parallel scheme, which can yield very high load balance and parallel efficiency. We have established a master-worker parallel scheme for the calculation of HFX forces as that for the construction of HFX matrix based on the dynamic parallel distribution algorithm (Shang et al., 2020). In this scheme, one message passing interface (MPI) process is designated as the master to manage the distribution of the shell quartets, while the remaining processes act as the workers responsible for integral evaluation. The computational task corresponding to the total amount of shell quartets is obtained by multiplying the size of shell-pair lists. Then, the shell quartets are assigned by the master to the worker processes in batches by request at a time. Each worker process requests individual and batched shell quartets from the master, and computes the ERI derivatives and HFX forces with integral screening. Once the current tasks are completed, the worker process continues to request new shell quartets until there are no tasks left. In order to further reduce the data communication for the density matrix, we replicate the global density matrix on each individual MPI process. Unlike the HFX matrix construction, the calculation of HFX forces does not require a MPI_ALLREDUCE operation, thus global communication can be avoided.

In recent years, heterogeneous architectures with dedicated accelerators have become increasingly available in modern HPC systems. As one of the most widely used accelerators, graphics processing unit (GPU) is designed specifically for handling massive parallelism and performing multiple computational tasks simultaneously. Nowadays, numerous quantum chemistry software packages have been equipped with GPU support, in particular, to accelerate the evaluation of ERIs and their subsequent contraction for constructing the HFX matrix (Ufimtsev and Martinez, 2008, 2009; Barca et al., 2020). The GPU acceleration can be accomplished by mapping ERIs onto GPU threads, using either a one-block-one-contracted-integral algorithm or a one-thread-one-contracted-integral algorithm (Ufimtsev and Martinez, 2008). Regarding the master-worker parallel distribution, we can set an appropriate batch size and map the batched ERIs of a request across the GPU blocks or threads, which depends on the available resources on the GPU associated with each worker process. Therefore, our parallelization strategy presented above could also be extended to CPU-GPU heterogeneous parallelism, even though such an extension is not covered in the present work.

In order to illustrate the parallel scalability of our implementation, we perform Γ-only HSE06 calculations for the Si crystal with a supercell containing 512 atoms by using different CPU cores. The default parameters are used for NAO2GTO fitting and integral screening. Our calculations are performed on Intel(R) Xeon(R) CPUs (6258R CPU@2.70GHz). The wall time for building the HFX matrix in the last SCF step is recorded, while the reported wall time for computing the HFX forces only includes the second term in Eq. 13

Figure 3A shows the change of wall times with respect to the number of CPU cores ranging from 48 to 768. The calculation of HFX forces takes 222.1 and 99.2 s for 48 and 768 CPU cores, respectively, which is 1.4–2.2 times faster than the HFX matrix construction (2142.9 and 1521.9 s). In the master-worker dynamic parallelization of HFX force calculation, the load balance can be effectively achieved, and only point-to-point communication of shell quartet indices is needed between the master and worker processes. Thus, the HFX force calculation scales nearly perfectly up to 768 CPU cores with a very high parallel efficiency of 95.9% as expected. However, the parallel efficiency for the construction of HFX matrix is significantly reduced to 57.8%. This reduction can be attributed to all-to-all communications required for building the global HFX matrix, which has been demonstrated in our previous work (Shang et al., 2020).

With such a good parallel scalability, we demonstrate the linear-scaling behavior of our implementation with respect to system size in parallel. We perform a series of Γ-only HSE06 calculations for the Si crystal with different supercells containing from 64 up to 1024 atoms on 240 CPU cores, in which the HFX matrix construction still maintains high parallel efficiency. As shown in Figure 3B, the wall time of both HFX matrix and force computations scale linearly with respect to system size. In particular, the linear-scaling calculation of HFX forces has a smaller prefactor than that of HFX matrix.

As a prototypical photocatalyst, TiO₂ is one of the most intensively studied materials, and polarons often play a decisive role in its applications (De Lile et al., 2022). For bulk rutile TiO₂, excess electrons can self-trap to form small polarons associated with local lattice distortion (Setvin et al., 2014). It has shown that hybrid functionals are sufficiently accurate to describe the formation and properties of small polarons in rutile TiO₂(Janotti et al., 2013; Elmaslmane et al., 2018; De Lile et al., 2022). Herein, we apply our code to investigate the small polaron due to the excess electron in bulk rutile TiO₂ with HSE06. In all our calculations, the experimental lattice constants of a = 4.594 Å and c = 2.959 Å are used. The calculated band gap for rutile TiO₂ is 3.28 eV, slightly higher than the experimental value of 3.03 eV (Amtout and Leonelli, 1995) but lower than the reported HSE06 value of 3.39 eV (Landmann et al., 2012). To simulate the formation of small polaron, a 3 × 3 × 4 supercell containing 216 atoms and −1|e| net charge is used for spin-polarized HSE06 calculations. The k-point meshes of 2 × 2 × 2 and 4 × 4 × 4 are chosen for structural relaxation and electronic structure calculations, respectively. One Ti atom is specified an initial displacement ( ∼ 0.18 Å) for localization of the polaron, and all atomic coordinates are relaxed until the forces acting on each atom are less than 0.04 eV/Å.

After full structural optimization, we observe a local lattice distortion around one Ti ion in the electron-doped rutile TiO₂. Compared to the pristine structure, the two Ti-O bonds perpendicular to the c-axis relax outward, increasing from 1.981 to 1.991 Å in length. In particular, the other four bonds with an initial length of 1.948 Å undergo two distinct changes: two of them increase to 2.011 Å, while the remaining two decrease to 1.891 Å. Figure 4A shows the band structure, where we can find a localized spin-electron state located at roughly 0.88 eV below the conduction band minimum (CBM). We also plot the spin density for this localized state in Figure 4B. As expected, the spin density is localized on the single Ti ion with a local lattice distortion, indicating the formation of a small electron polaron. Our results are in good agreement with the reported HSE06 results, in which an electron polaron state at 0.77 eV below the CBM was predicted by using VASP (Janotti et al., 2013).