Hybrid functionals for periodic systems in the density functional tight-binding method

TL;DR

This paper extends the density functional tight-binding (DFTB) method to include screened range-separated hybrid functionals, enabling accurate and efficient modeling of periodic systems with improved electronic property predictions.

Contribution

It generalizes the theoretical foundation of DFTB to incorporate SRSH functionals with periodic boundary conditions, including techniques for handling Fock exchange in reciprocal space.

Findings

Accurate prediction of polarization-induced gap renormalization in molecular crystals.

Demonstrated convergence and efficiency for polyacenes and layered materials.

Reduced computational cost compared to first principles methods.

Abstract

Screened range-separated hybrid (SRSH) functionals within generalized Kohn-Sham density functional theory (GKS-DFT) have been shown to restore a general $1/(r\varepsilon)$ asymptotic decay of the electrostatic interaction in dielectric environments. Major achievements of SRSH include an improved description of optical properties of solids and correct prediction of polarization-induced fundamental gap renormalization in molecular crystals. The density functional tight-binding method (DFTB) is an approximate DFT that bridges the gap between first principles methods and empirical electronic structure schemes. While purely long-range corrected RSH are already accessible within DFTB for molecular systems, this work generalizes the theoretical foundation to also include screened range-separated hybrids, with conventional pure hybrid functionals as a special case. The presented formulation and implementation is also valid for periodic boundary conditions (PBC) beyond the $\Gamma$-point. To treat periodic Fock exchange and its integrable singularity in reciprocal space, we resort to techniques successfully employed by DFT, in particular a truncated Coulomb operator and the minimum image convention. Starting from the first principles Hartree-Fock operator, we derive suitable expressions for the DFTB method, using standard integral approximations and their efficient implementation in the DFTB+ software package. Convergence behavior is investigated and demonstrated for the polyacene series as well as two- and three-dimensional materials. Benzene and pentacene molecular and crystalline systems show the correct polarization-induced gap renormalization by SRSH-DFTB at heavily reduced computational cost compared to first principles methods.

Equations

$$\frac{1}{r}=\underset{\text{DFT}}{\underbrace{\frac{1-(\alpha +\beta )+\beta {\mathrm{e}}^{-\omega r}}{r}}}+\underset{\text{HF}}{\underbrace{\frac{\alpha +\beta (1-{\mathrm{e}}^{-\omega r})}{r}}},$$

$$E[\rho ]$$

$$+[1-(\alpha +\beta )]E_{\mathrm{x}}^{\mathrm{DFT}}+\alpha E_{\mathrm{x}}^{\mathrm{HF}}+\beta (E_{\mathrm{x},\mathrm{sr}}^{\omega ,\mathrm{DFT}}+E_{\mathrm{x},\mathrm{lr}}^{\omega ,\mathrm{HF}}),$$

$$E_{\mathrm{xc},\mathrm{loc}}[\rho ]$$

$$E$$

$$+\frac{N}{2}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}}{w}_{{\mathbf{k}}^{\prime }}\sum _{ij}{f}_{i\mathbf{k}}{f}_{j{\mathbf{k}}^{\prime }}\iint \frac{|{\psi }_{i\mathbf{k}}(\mathbf{r}){|}^2|{\psi }_{j{\mathbf{k}}^{\prime }}({\mathbf{r}}^{\prime }){|}^2}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }$$

$$-\frac{N}{4}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}}{w}_{{\mathbf{k}}^{\prime }}\sum _{ij}{f}_{i\mathbf{k}}{f}_{j{\mathbf{k}}^{\prime }}\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }$$

$$\times {\psi }_{i\mathbf{k}}^{\ast }(\mathbf{r}){\psi }_{j{\mathbf{k}}^{\prime }}(\mathbf{r})\frac{\alpha +\beta \left[1-{\mathrm{e}}^{-\omega |\mathbf{r}-{\mathbf{r}}^{\prime }|}\right]}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\psi }_{j{\mathbf{k}}^{\prime }}^{\ast }({\mathbf{r}}^{\prime }){\psi }_{i\mathbf{k}}({\mathbf{r}}^{\prime })$$

$$+E_{\mathrm{xc},\mathrm{loc}}[\rho ]+E_{\mathrm{NN}}.$$

$${\beta }_{\mu }^{\mathbf{k}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum _{\mathbf{g}}{\varphi }_{\mu }(\mathbf{r}-\mathbf{g}){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{g}}.$$

$${\psi }_{i\mathbf{k}}(\mathbf{r})$$

$$O_{\mu \nu }(\mathbf{g})$$

$$O_{\mu \nu }(\mathbf{k})$$

$$P_{\mu \nu }(\mathbf{k})$$

$$P_{\mu \nu }(\mathbf{g})$$

$$\left({\varphi }_{\mu }^{\mathbf{l}}{\varphi }_{\nu }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\lambda }^{\mathbf{g}}\right):=\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }$$

$$\times {\varphi }_{\mu }^{\ast }(\mathbf{r}-\mathbf{l}){\varphi }_{\nu }(\mathbf{r}-\mathbf{h})\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\varphi }_{\kappa }^{\ast }({\mathbf{r}}^{\prime }-\mathbf{m}){\varphi }_{\lambda }({\mathbf{r}}^{\prime }-\mathbf{g})$$

$${\left({\varphi }_{\mu }^{\mathbf{l}}{\varphi }_{\nu }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\lambda }^{\mathbf{g}}\right)}^{\mathrm{lr},\omega }:=\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }$$

$$\times {\varphi }_{\mu }^{\ast }(\mathbf{r}-\mathbf{l}){\varphi }_{\nu }(\mathbf{r}-\mathbf{h})\frac{1-{\mathrm{e}}^{-\omega |\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\varphi }_{\kappa }^{\ast }({\mathbf{r}}^{\prime }-\mathbf{m}){\varphi }_{\lambda }({\mathbf{r}}^{\prime }-\mathbf{g}),$$

$$E$$

$$+\frac{1}{2}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}}{w}_{{\mathbf{k}}^{\prime }}\sum _{\mu \nu \kappa \lambda }{P}_{\nu \mu }(\mathbf{k})P_{\lambda \kappa }({\mathbf{k}}^{\prime })$$

$$\times \sum _{\mathbf{g}\mathbf{h}\mathbf{m}}{\mathrm{e}}^{i{\mathbf{k}}^{\prime }(\mathbf{g}-\mathbf{m})}\left({\varphi }_{\mu }^0{\varphi }_{\nu }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\lambda }^{\mathbf{g}}\right){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{h}}$$

$$-\frac{1}{4}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}}{w}_{{\mathbf{k}}^{\prime }}\sum _{\mu \nu \kappa \lambda }{P}_{\nu \mu }(\mathbf{k})P_{\lambda \kappa }({\mathbf{k}}^{\prime })\sum _{\mathbf{g}\mathbf{h}\mathbf{m}}{\mathrm{e}}^{-i{\mathbf{k}}^{\prime }(\mathbf{m}-\mathbf{h})}$$

$$\times \left[\alpha \left({\varphi }_{\mu }^0{\varphi }_{\lambda }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\nu }^{\mathbf{g}}\right)+\beta {\left({\varphi }_{\mu }^0{\varphi }_{\lambda }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\nu }^{\mathbf{g}}\right)}^{\mathrm{lr},\omega }\right]{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{g}}.$$

$$\sum _{\nu }{c}_{\nu i}(\mathbf{k})\left[H_{\mu \nu }(\mathbf{k})-{\epsilon }_i{S}_{\mu \nu }(\mathbf{k})\right]=0,\forall i$$

$$H_{\mu \nu }(\mathbf{k})$$

$$+\sum _{{\mathbf{k}}^{\prime }}{w}_{{\mathbf{k}}^{\prime }}\sum _{\lambda \kappa }{P}_{\lambda \kappa }({\mathbf{k}}^{\prime })\sum _{\mathbf{g}\mathbf{h}\mathbf{m}}{\mathrm{e}}^{i{\mathbf{k}}^{\prime }(\mathbf{g}-\mathbf{m})}$$

$$\times \left({\varphi }_{\mu }^0{\varphi }_{\nu }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\lambda }^{\mathbf{g}}\right){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{h}}$$

$$-\frac{1}{2}\sum _{{\mathbf{k}}^{\prime }}{w}_{{\mathbf{k}}^{\prime }}\sum _{\lambda \kappa }{P}_{\lambda \kappa }({\mathbf{k}}^{\prime })\sum _{\mathbf{g}\mathbf{h}\mathbf{m}}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{g}}{\mathrm{e}}^{i{\mathbf{k}}^{\prime }(\mathbf{h}-\mathbf{m})}$$

$$\times \left[\alpha \left({\varphi }_{\mu }^0{\varphi }_{\lambda }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\nu }^{\mathbf{g}}\right)+\beta {\left({\varphi }_{\mu }^0{\varphi }_{\lambda }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\nu }^{\mathbf{g}}\right)}^{\mathrm{lr},\omega }\right],$$

$$v_{\mu \nu }^{\mathrm{xc},\mathrm{loc}}[\rho ](\mathbf{g})=\int {\varphi }_{\mu }^{\mathbf{g}}(\mathbf{r}{)}^{\ast }{v}_{\mathrm{xc},\mathrm{loc}}[\rho ](\mathbf{r}){\varphi }_{\nu }^0(\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$E_{\mathrm{xc},\mathrm{loc}}[{\rho }_0+\delta \rho ]=E_{\mathrm{xc},\mathrm{loc}}[{\rho }_0]+\frac{1}{N}\int v_{\mathrm{xc},\mathrm{loc}}[{\rho }_0](\mathbf{r})\delta \rho (\mathbf{r})\mathrm{d}\mathbf{r}$$

$$+\frac{1}{2N}\iint f_{\mathrm{xc},\mathrm{loc}}[{\rho }_0](\mathbf{r},{\mathbf{r}}^{\prime })\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }+\mathcal{O}(\delta {\rho }^3),$$

$$E_{\mathrm{xc},\mathrm{loc}}[\rho ]$$

$$+\frac{1}{2}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}}{w}_{{\mathbf{k}}^{\prime }}\sum _{\mu \nu \kappa \lambda }\Delta P_{\nu \mu }(\mathbf{k})\Delta P_{\lambda \kappa }({\mathbf{k}}^{\prime })$$

$$\times \sum _{{\mathbf{g}}^{\prime }}\sum _{\mathbf{l}{\mathbf{l}}^{\prime }}{\mathrm{e}}^{i\mathbf{k}\cdot {\mathbf{g}}^{\prime }}{\mathrm{e}}^{i{\mathbf{k}}^{\prime }({\mathbf{l}}^{\prime }-\mathbf{l})}{f}_{\mu \nu \kappa \lambda }^{\mathrm{xc},\mathrm{loc}}(0,{\mathbf{g}}^{\prime },\mathbf{l},{\mathbf{l}}^{\prime })$$

$$f_{\mu \nu \kappa \lambda }^{\mathrm{xc},\mathrm{loc}}(\mathbf{g},{\mathbf{g}}^{\prime },\mathbf{l},{\mathbf{l}}^{\prime })$$

$$\times {\varphi }_{\mu }^{\mathbf{g}}(\mathbf{r}{)}^{\ast }{\varphi }_{\nu }^{{\mathbf{g}}^{\prime }}(\mathbf{r})f_{\mathrm{xc},\mathrm{loc}}(\mathbf{r},{\mathbf{r}}^{\prime }){\varphi }_{\kappa }^{\mathbf{l}}({\mathbf{r}}^{\prime }{)}^{\ast }{\varphi }_{\lambda }^{{\mathbf{l}}^{\prime }}({\mathbf{r}}^{\prime }),$$

$$E$$

$$E_{\text{DFT}}^{(2)}$$

$$\times {\mathrm{e}}^{i{\mathbf{k}}^{\prime }(\mathbf{g}-\mathbf{m})}\left[\left({\varphi }_{\mu }^0{\varphi }_{\nu }^{\mathbf{h}}|{\varphi }_{\kappa }^{\mathbf{m}}{\varphi }_{\lambda }^{\mathbf{g}}\right)+f_{\mu \nu \kappa \lambda }^{\mathrm{xc},\mathrm{loc}}(0,\mathbf{h},\mathbf{m},\mathbf{g})\right]$$

$$E_{\text{HF}}^{(2)}$$

Full text

I Introduction

As a semi-empirical method, density functional tight binding (DFTB) [1, 2] fills the gap between methods such as Hartree-Fock [3, 4] or Kohn-Sham density functional theory (DFT) [5, 6] and fully empirical force-fields in the domain of computational chemistry, condensed matter physics and materials science. Its high ratio of accuracy to computational cost renders DFTB well-suited for extended systems with large unit cells or long timescale molecular dynamics (MD). Over the last two decades, the original DFTB formalism by Seifert et al. [1] has been expanded by a number of extensions, including self-consistent charge SCC-DFTB [2] and its extension to third order DFTB3 [7], spin and spin-orbit interactions [8], time dependent TD-DFTB [9], real-time rTD-DFTB using propagation of the reduced one body density matrix and Ehrenfest dynamics [10, 11], machine learning enhanced repulsive potentials [12, 13] as well as non-equilibrium Green’s function based electron transport [14].

For molecular calculations, range-separated hybrid functionals (RSH) [15, 16, 17, 18, 19] which constitute a mixture of non-local Fock like and (semi-)local exchange of DFT have been established as standard technique to address the inherent electronic self-interaction error (SIE) [20] of DFT and restore the piecewise-linear [21, 22] behavior of the exact exchange-correlation functional between integer occupations. Niehaus and Della Sala [23] generalized the DFTB formalism to non-periodic long-range corrected hybrid functionals (LC-DFTB), based on GKS-DFT and the density matrix as basic variable in the expansion of the Kohn-Sham energy functional, which was later implemented in the DFTB+ [24, 25] software package.

For periodically repeating structures, difficulties due to the Coulomb singularity of Fock exchange initially prevented an immediate and widespread adoption of hybrid functionals for solids. The CRYSTAL [26] software package, based on the work of Pisani and Dovesi [27], provided the first publicly available implementation of periodic Fock exchange, paving the way for making periodic hybrid functionals readily accessible for solids. In recent years effort has been made to develop reliable schemes for treating the singularity, including pioneering work by Gygi and Baldereschi [28] who lifted the singularity by introducing auxiliary functions, Spencer and Alavi [29] resorting to a truncated Coulomb operator that is relatively simple to implement and does not posses a singularity in reciprocal space, and most recently Sundararaman and Arias’ [30] analytical proof of Wigner-Seitz truncation as an ideal method for regularizing the Coulomb potential in the exchange kernel.

General range-separated hybrid functional implementations have been developed for plane-wave [31, 32, 33, 34, 35] and localized numerical [36, 37, 38] or Gaussian-type [39, 40, 41, 42] orbital based DFT codes. Here we derive suitable expressions for the DFTB method and implement a real-space formulation of periodic Fock exchange in the DFTB+ software package. We compare two algorithmic solutions that follow the footsteps of the truncated Coulomb interaction (TCI) by Guidon et al. [41] as well as the minimum image convention (MIC) by Tymczak et al. [42] that utilizes the full, unscreened, Coulomb interaction.

The present work is structured as follows. In Section II, we outline the basic periodic hybrid-functional DFTB theory, starting from the expansion of the total energy functional, deriving the periodic Fock exchange Hamiltonian within DFTB (by imposing Born–von Kármán (BvK) periodic boundary conditions for the density matrix of a finite mesh of k-points). Section III translates the equations into a formulation that is compatible with the neighbor list concept of DFTB+, and additionally gives expressions optimized for $\mathrm{\Gamma}$-point only k-point sampling, enabling efficient simulation of large systems with thousands of atoms. This Section also contains benchmarks covering the scaling of the $\mathrm{\Gamma}$-point implementation with system size and the parallel performance of the k-point implementation. Total energy and band-gap convergence behavior is then investigated and demonstrated by Section IV for the polyacene series, complemented by armchair graphene nanoribbons, two-dimensional h-BN monolayers and GaAs bulk in the Supplemental Material [43]. In Section V, benzene and pentacene molecular and crystalline systems are shown to exhibit the correct polarization-induced gap renormalization, that occurs if the surrounding dielectric medium is properly taken into account. We close this work with a summary of our findings and by providing a brief outlook in Section VI.

We should note that the focus of the present work lies in the theory and implementation of the method, rather than an in-depth benchmarking of its accuracy. The latter will be the topic of a forthcoming article.

II Theory

II.1 Periodic GKS-DFT formalism

In the following we outline the basic periodic KS-DFT formalism, including Fock exchange. This provides the notation and forms the basis for the DFTB approach described in Section II.2. All quantities are given in atomic units throughout and we will denote the crystal momentum as $\bm{k}$ and $\bm{k}^{\prime}$, while real space lattice vectors are denoted as $\bm{g}$, $\bm{h}$, $\bm{l}$ and $\bm{m}$.

According to the Coulomb-attenuating method (CAM) [18], based on pioneering works by Gill [44] and Savin [45], the electron-electron interaction is partitioned into short- and long-range contributions using the adiabatic connection theorem [46]

[TABLE]

where the short-range part is handled by a modified purely density functional, that ensures mutual error cancellation of local exchange and correlation, while the exact Fock exchange enforces the correct asymptotic decay in the long-range limit. The parameters $\alpha,\beta$ and $\omega$ determine the fraction of global and long-range exact Fock exchange, as well as the value of the smooth range-separation function, which we have assumed to be of Yukawa type.

Complementing the usual kinetic energy $T_{0}[\rho]$ of the auxiliary system of non-interacting electrons in GKS-DFT with the classical Coulomb interaction $E_{\mathrm{H}}[\rho]$ (Hartree term), external potential $v^{\mathrm{ext}}(\bm{r})$ and nuclear-repulsion energy $E_{\mathrm{NN}}$, the above partitioning leads to the total energy expression per unit cell (UC)

[TABLE]

with $N$ the total number of such cells in the crystal. This expression can be shortened by introducing the local part of the xc-functional

[TABLE]

Note that $E_{\mathrm{xc,loc}}[\rho]=E_{\mathrm{xc,loc}}^{\alpha,\beta,\omega}[\rho]$, but we refrain from explicitly stating the additional parameter dependencies for brevity. Expressing the individual contributions of Eq. (II.1) in terms of orbitals $\psi_{i\bm{k}}(\bm{r})$, with additional quantum number $\bm{k}$, the total energy reads

[TABLE]

We describe a spin-unpolarized formalism for closed shell systems, where the spin degrees of freedom have already been summed up. The occupation of eigenstate $i$ at crystal momentum $\bm{k}$ is denoted as $f_{i\bm{k}}\in[0,2]$. The number of states $i$ is given by the number of basis functions in the UC. We further introduce weights $w_{\bm{k}}$, with normalization $\sum_{\bm{k}}w_{\bm{k}}=1$, that arise from sampling the Brillouin zone by only selecting a subset of all wave-vectors compatible with the BvK cell.

For a crystal that is invariant with respect to cell translations, the GKS-orbitals are expected to obey the Bloch theorem and extend throughout the whole crystal. We introduce Bloch-functions $\beta_{\mu}^{\bm{k}}(\bm{r})$, that emerge from a unitary transformation of the atomic orbitals $\phi_{\mu}(\bm{r})$ for orbital $\mu$ centered on an atom in the reference cell, when shifted by any real-space lattice vector $\bm{g}$:

[TABLE]

Expressing the wavefunctions as a linear combination of crystalline orbitals leads to states that equally possess Bloch-wave character, therefore satisfying the Bloch-condition $\psi_{i\bm{k}}(\bm{r}+\bm{g})=\psi_{i\bm{k}}(\bm{r})\mathrm{e}^{i\bm{k}\cdot\bm{g}}$

[TABLE]

with eigenvector coefficients $c_{\mu i}(\bm{k})$ attributed to orbital $\mu$ and eigenstate $i$. By exploiting the translational symmetry, $\hat{\mathrm{O}}(\bm{r}+\bm{g},\bm{r}^{\prime}+\bm{g})=\hat{\mathrm{O}}(\bm{r},\bm{r}^{\prime})$, of an operator $\hat{\mathrm{O}}$, folding operations for transitions from direct to reciprocal space (and vice versa) are obtained:

[TABLE]

The special case $\hat{\mathrm{O}}=\mathbf{1}$ refers to the overlap matrix elements $S_{\mu\nu}$. We introduce the convention that the real-space shifts in the arguments of Hamiltonian and overlap refer to the first orbital, while the second remains in the reference cell. In reciprocal space, the density matrix is built from the eigenvector coefficients and occupations

[TABLE]

while transformations according to Eq. (7) and Eq. (8) apply.

For the two-electron, four-center integrals we resort to Mulliken’s notation and distinguish between an unscreened and screened, long-range Coulomb kernel

[TABLE]

with shorthand notation $\phi_{\mu}^{\bm{g}}:=\phi_{\mu}(\bm{r}-\bm{g})$. We represent the total energy of Eq. (II.1) in terms of the new Bloch basis and exploit the translational symmetry of an arbitrary (in general, non-local) operator in direct space to choose the reference (zeroth) cell and carry out a trivial lattice summation, resulting in

[TABLE]

Here, $h_{\mu\nu}(\bm{k})$ refers to the one-electron integral matrix elements of the kinetic energy and external potential.

By applying the variational principle to Eq. (II.1) we obtain the secular equation

[TABLE]

with $\bm{k}$-dependent Hamiltonian matrix elements

[TABLE]

where the local part of the exchange-correlation potential (i.e., the functional derivative of $E_{\mathrm{xc,loc}}$) is introduced as

[TABLE]

Bloch-functions are bases for the irreducible representations of the translation group and super-matrices in reciprocal space, such as $H_{\mu\nu}(\bm{k})$, can be transformed into a block-diagonal form by unitary transformation [47]. This gives rise to an independent diagonalization of Eq. (15) for each different k-point.

II.2 Density functional tight binding

In the following we outline the periodic RSH-DFTB formalism of second order in the energy expansion.

The minimal valence-only basis set $\{\phi_{\mu}\}$ of atomic orbitals entering the eigenfunction ansatz of Eq. (6) is obtained by performing first principles RSH calculations for neutral and spin-unpolarized pseudo-atoms, as described in Refs. 7, 24. We continue by approximating the local part of the exchange-correlation functional, by expanding around the reference up to second order in the perturbation

[TABLE]

thus linearizing the exchange-correlation potential. The density matrix of Eq. (9) is decomposed into reference and perturbation, such that ${\bf P}={\bf P}^{(0)}+\Delta{\bf P}$. Usually, the reference ${\bf P}^{(0)}=\sum_{A}{\bf P}_{A}$ is constructed as a superposition of densities of non-interacting atoms, where the sum runs over all atoms ($A$) in the unit cell. Representing Eq. (18) in the Bloch basis of Eq. (6) then yields

[TABLE]

where the matrix elements of first and second functional derivatives are denoted as $v_{\mu\nu}^{\mathrm{xc,loc}}$ and $f_{\mu\nu\kappa\lambda}^{\mathrm{xc,loc}}$, respectively.

We continue by inserting Eq. (19) into the original energy functional of Eq. (II.1), leading to

[TABLE]

with $E_{\mathrm{rep}}$ covering all terms that solely depend on the reference density and $E^{(2)}_{\text{DFT}}$ and $E^{(2)}_{\text{HF}}$ being of second order in $\Delta{\bf P}$

[TABLE]

Further, Eq. (21) introduced the zeroth-order Hamiltonian $H^{(0)}_{\mu\nu}(\bm{k})$, that is defined as follows

[TABLE]

As usual in the DFTB framework, we adopt the two-center approximation and replace onsite-blocks in the Hamiltonian by diagonal matrices with free atom eigen-energies, $\epsilon_{\mu}^{\text{free}}$, ensuring the correct limit on dissociation and leading to

[TABLE]

$\rho_{A}=\rho_{A}(\bm{r}-\bm{g})$ and $\rho_{B}=\rho_{B}(\bm{r})$ denote atomic densities, derived from appropriate pseudo-atom calculations. $\epsilon_{\mu}^{\text{free}}$, as well as the non-diagonal elements of Eq. (26), are obtained from RSH-DFT calculations and latter stored for high-symmetry orbital configurations as a function of distance between atoms $A$ and $B$ in Slater-Koster tables [48]. The matrix elements $H_{\mu\nu}^{(0)}(\bm{k})$ are recovered by transforming according to Eq. (8).

In line with conventional DFTB and as generalized to periodic boundary conditions, the repulsive energy is approximated by a sum of fast decaying pair potentials [49]

[TABLE]

either determined by a higher level of theory [7] or fitted to empirical data [50]. The sums over atoms $A,B$ are restricted to the reference unit cell, while an additional direct lattice sum $\bm{g}$ also accounts for contributions of images in neighboring cells.

The term $E^{(2)}_{\text{DFT}}$, defined in Eq. (22) is treated by applying the Mulliken approximation

[TABLE]

to the four-center integrals

[TABLE]

In line with the convention introduced by Eq. (7), the real-space shifts in the arguments of $\gamma$ are associated with the first orbital, while the second orbital remains in the central cell. We follow the reasoning of conventional second-order DFTB, by approximating the orbital products $|\phi_{\mu}|^{2}$ as exponentially decaying spherically symmetric charge densities, leading to three integral parameterizations

[TABLE]

with analytical expressions [2, 24]. A distinction is made between screened long-range (lr) and unscreened full-range (fr) kernels. The parameter $\tau_{A}$ is obtained from requiring RSH-DFT and RSH-DFTB yield the same second derivative of the total energy for an atom with respect to orbital occupation, i.e., predicting the same chemical hardness [24]. A detailed derivation is provided in Appendix A.

Now we have all the prerequisites to treat the semi-local energy contribution $E^{(2)}_{\text{DFT}}$. To this end we compute Mulliken populations, $q_{\mu}$, of orbital $\mu$

[TABLE]

evaluated in direct or reciprocal space. The net charges of atom $A$ are obtained by comparing the populations with the neutral atom ($Z_{A}$)

[TABLE]

The final result is obtained by summing over all atoms $A,B$ in the unit cell and accounting for any periodic images of $B$, as captured by the sum over direct lattice vectors $\bm{g}$ extending throughout the crystal

[TABLE]

II.3 Periodic Fock exchange in DFTB

To complete the theoretical foundation of periodic RSH-DFTB, the energy contributions due to the additional Fock terms require special treatment. We restrict the derivation to the screened, long-range Coulomb kernel, since its full-range counterpart is contained as the limiting case of a large value of the range-separation parameter, $\omega$.

Firstly, the (back) Fourier transformation of Eq. (11) is identified in order to simplify the expression, leading to

[TABLE]

We now apply the Mulliken approximation and arrive at an expression that is compatible with the DFTB formalism

[TABLE]

For later implementation using a neighbor list, as discussed in Section III.1, it is advantageous to carry out an index shift $\bm{l}\longrightarrow\bm{g}-\bm{h}-\bm{l}$ that simplifies the arguments of the real-space overlaps

[TABLE]

By applying the variational principle and writing the above expression with respect to the density matrices, the corresponding ground-state Hamiltonian emerges as

[TABLE]

This illustrates that the evaluation of the associated energy contribution becomes straightforward once the exchange Hamiltonian $\Delta H^{x,\mathrm{lr}}_{\mu\nu}(\bm{k})$ is available

[TABLE]

Deriving the energy and Hamiltonian for global Fock-like exchange, instead of a range-separated expression, simply requires substituting $\gamma_{\mu\nu}^{\mathrm{lr,HF}}$ with $\gamma_{\mu\nu}^{\mathrm{fr,HF}}$.

III Implementation

The extension of the DFTB formalism, according to Sections II.2 and II.3, requires modifications to a) the parameterization suite skprogs [51] and b) the main DFTB+ [25] code. The newly developed routines are publicly available in a) the main branch of the official repository and b) a pull request of the respective development branch to the main branch of the official repository [52], currently under code review by the maintainers.

For the zeroth-order Hamiltonian construction, we generalized the scheme of Lutsker and co-workers [24] to handle global Hartree-Fock exchange in addition to the already available screened kernels. This enables pre-tabulation of the $H^{(0)}$ and $S$ matrix elements for general CAM xc-functionals and further includes a generalized scheme to determine the decay constants according to Eq. (61).

The central performance critical task of developing an efficient implementation to construct the Fock exchange matrix, as provided by Eq. (II.3), is based on the neighbor-list-based design of DFTB+. This utilizes the sparsity [53] pattern induced by the spatial decay of the real-space overlap matrix elements in the Slater-Koster tables.

III.1 Neighbor list based algorithm

In order to re-formulate Eq. (II.3), utilizing the concept of neighbor lists, we first introduce the notation $\bar{L}\in\mathcal{N}(M)$ which refers to atom $\bar{L}$ being a neighbor of atom $M$, where atom $M$ is located in the central cell but $\bar{L}$ could be inside a periodically repeated neighboring unit cell, with $L$ being the corresponding atom in the central cell. In a similar fashion, the atomic orbitals $\bar{\lambda}[\bar{L}]$ are not restricted to atoms in the central cell, but the corresponding orbital in the central cell is labeled as $\lambda[L]$. This enables us to drop the cell index of the real-space overlap matrices, since they are implicitly included for atoms located outside the central cell.

We also follow DFTB+ specific conventions, such that the real-space shift of the overlap $S$ refers to its first index, in line with the notation of Section II.1, and in particular Eq. (7). The same reasoning also applies to the real-space density matrix elements. Let $\lambda$ denote an orbital that is folded back into the central cell from a periodic image, we may then write $S_{\bar{\lambda}\mu}=S_{\lambda\mu}(\bm{h})$, leading to

[TABLE]

Eq. (III.1) exploits the fact that, except in the case of shell-resolved DFTB, the $\gamma$-functions depend only on the atomic species and not the individual orbitals, i.e. $\gamma^{\mathrm{lr,HF}}_{\bar{\lambda}\bar{\kappa}}(\bm{g})=\gamma^{\mathrm{lr,HF}}_{\bar{L}\bar{K}}(\bm{g})$ with $\bar{\lambda}[\bar{L}]$ and $\bar{\kappa}[\bar{K}]$. Analytic expressions for atomic forces, i.e. the negative derivative of Eq. (III.1) with respect to the ion positions, have been derived and implemented as well. Although the force expressions for the general $\bm{k}$-point implementation have already been validated against numerical derivatives, their algorithmic optimization is subject to ongoing development and the current implementation is of limited use for production applications. However, optimized forms of Eq. (III.1) for calculations restricted to the $\Gamma$-point, including the energy gradients, are provided in Section B.

So far, two of the three previously infinite lattice summations have been replaced by well-defined finite summations over the neighbor list, leaving the yet unbounded $\bm{g}$-summation to discuss.

Practical calculations employ a finite set of k-points to sample the first Brillouin zone, with $N_{k}$ unit cells spanning the BvK supercell. Restricting the Bloch basis to a finite BvK supercell leads to finite-size errors, since it is not complete with respect to all possible wavevectors, $\bm{k}$, of the infinite crystal. The density matrix as introduced in Eq. (11) is by construction BvK periodic: $P_{\mu\nu}(\bm{g})=P_{\mu\nu}(\bm{g}+\bm{G})$, with $\bm{G}$ denoting a BvK super-lattice vector. This, means it repeats at the boundaries of the BvK supercell without a phase factor. As shown by Irmler et al. [39], this artificial periodicity causes lattice sums in the Fock exchange energy expression to diverge. We have implemented two of the widespread schemes to remedy this issue, namely a truncated Coulomb interaction [29] (TCI) and an adaption of the minimum image convention (MIC) [42]. These enable a robust implementation of periodic Fock exchange. A more pleasant consequence of the BvK periodicity is that the real-space formulation only requires $N_{k}$ density matrices, $P_{\mu\nu}(\bm{g})$, to be kept in storage. This is because every lattice shift, $\bm{g}$, may be folded back into the central BvK cell.

III.2 Truncated Coulomb interaction

In RSH-DFTB, truncating the Coulomb kernel of the four-center integrals is equivalent to limiting the range of the $\gamma$-integrals

[TABLE]

with an adjustable real-space cutoff radius $R_{\mathrm{c}}$. To avoid interactions with the neighboring BvK supercells in simple-cubic systems, Spencer and Alavi [29] linked $R_{\mathrm{c}}$ to the number of k-points (which determine the BvK supercell volume). A more robust scheme for arbitrary lattice geometries, as implemented in this work, is to determine the maximum radius of a sphere that still fits within the BvK supercell [39].

During self-consistent cycles the geometry and therefore also the $\gamma$-integrals do not change. Our implementation pre-tabulates all non-vanishing $\gamma_{MN}(R_{MN})$ within the cutoff sphere to speed up the Hamiltonian construction.

III.3 Minimum image convention

Another way of preventing divergent lattice sums was suggested by Tymczak et al. [42] by restricting the sum over super-lattice vectors according to the minimum image convention. Later, Irmler and co-workers [39] generalized this to arbitrary k-points. In this scheme the Coulomb interaction is unaltered and fully taken into account.

We adopt this idea by restricting the argument of the density matrix $\Delta P_{\mu\nu}(\bm{g})$ of Eq. (III.1), such that it does not involve orbitals outside its Wigner-Seitz cell. This naturally restricts the $\bm{g}$-summation which, depending on the size of the BvK cell, in turn depends on the k-point sampling employed. To determine the unit cells within the Wigner-Seitz cell of the BvK cell we employ an algorithm that does not assume a specific lattice geometry and works for arbitrary (linearly independent) lattice vectors.

III.4 Integral pre-screening

Regardless of whether the TCI or MIC algorithms are used, integral pre-screening targeting the density and overlap matrices has the potential to drastically reduce the computational cost of constructing the exchange Hamiltonian.

In direct self-consistent-field [54, 55] calculations, the Hamiltonian is often constructed iteratively. Following Lutsker and co-workers [24], the linearity of the Hamiltonian with respect to the density matrix allows representation of the Hamiltonian at the $n$-th self-consistent iteration as a sum of the Hamiltonian at the previous iteration $\Delta{\bf H}(\Delta{\bf P}^{n-1})$ plus a correction $\Delta{\bf H}(\Delta_{n}(\Delta{\bf P}))$. The change in the density matrix with respect to the previous iteration is denoted as $\Delta_{n}(\Delta{\bf P})$, resulting in

[TABLE]

This approach therefore exploits the rapid decay of $\Delta_{n}(\Delta{\bf P})$ with increasing cycles of self-consistency. During the Hamiltonian construction, matrix-matrix products of the form

[TABLE]

occur. The upper bound of Eq. (47) is provided by taking the individual absolute values of the factors

[TABLE]

where maximum estimates for the overlap $S^{\mathrm{max}}_{\bar{L}M}:=\max_{\bar{\lambda}[\bar{L}],\mu[M]}(|S_{\bar{\lambda}\mu}|)$ and density matrix $\Delta_{n}(\Delta{\bf P})^{\mathrm{max}}:=\max_{\bm{g}}(|\Delta_{n}(\Delta{\bf P}(\bm{g}))|)$ have been defined. If $S^{\mathrm{max}}_{\bar{L}M}S^{\mathrm{max}}_{\bar{K}N}\Delta_{n}(\Delta{\bf P})^{\mathrm{max}}<\varepsilon_{\mathrm{screen}}$, where the trivial summations over orbitals $\bar{\lambda},\bar{\kappa}$ have been absorbed by the integral screening parameter $\varepsilon_{\mathrm{screen}}$, the evaluation of the corresponding diatomic sub-block of $\Delta{\bf H}(\Delta_{n}(\Delta{\bf P}))$ is omitted. In preparation for the evaluation of Eq. (III.1), all occurring ${\bf S}(\Delta{\bf P}){\bf S}$-products are estimated and the terms requiring an explicit evaluation are distributed to available processors (provided that a message passing interface (MPI) parallelized version of DFTB+ is being used).

III.5 Scaling with system size

While supercells of several hundreds of atoms are often unattainable for proper long-range corrected hybrid functionals within RSH-DFT, we demonstrate that such cases are well within reach of RSH-DFTB, even on a single processor core and for relatively densely packed materials such as GaAs. Figure 1 compares the total wall-clock time of a RSH-DFTB $\Gamma$-point calculation using the LCY-PBE functional and Yukawa type range-separation function, to a conventional PBE-parameterized DFTB (referred to as PBE-DFTB) run. GaAs is computationally challenging due to its large number of interacting neighbors. The computational cost of LCY-PBE-DFTB turns out to be considerably higher than for conventional PBE-DFTB. However, considering that the benchmark was performed on a single CPU core only, calculations of large supercells with roughly 1000 atoms can be accomplished in reasonable time.

The higher cost of LCY-PBE-DFTB mainly originates from three facts: a) the pre-tabulation of all $\tilde{\gamma}_{MN}^{\mathrm{TC}}$ to build the super-matrix $\tilde{\gamma}^{\mathrm{TC}}$ (see Section B), b) the actual time spend on constructing the Hamiltonian of Eq. (B) and c) a higher number of total self-consistency steps, compared to traditional DFTB. Reason c) is expected, since (semi-)local DFTB requires only self-consistency with respect to the Mulliken populations, while RSH-DFTB introduces terms that depend on the full density matrix. Its self-consistency is with respect to the (real-space) matrices $\Delta P_{\mu\nu}(\bm{g})$. In other words, traditional DFTB mixes the in- and out-put Mulliken populations to propagate the self-consistent cycles, while RSH-DFTB mixes $\Delta P_{\mu\nu}(\bm{g})$, which proves to be more challenging and leads in most cases to an increased number of self-consistency steps, and therefore diagonalizations of the total Hamiltonian.

III.6 Parallel performance

The relatively high cost of constructing the exchange Hamiltonian in RSH-DFTB requires an efficient parallelization of this step, in order to exploit modern computing infrastructures and HPC facilities. While Section III.5 already demonstrated the suitability for large supercells, another common task is the calculation of smaller systems with a dense k-point sampling. Figure 2 illustrates the parallel performance for the energy evaluation of a primitive GaAs unit cell.

The scaling of the total wall-clock time with MPI processes is quite satisfactory and only saturates from about 100 cores onwards, for this system. One bottleneck that causes the parallel efficiency to drop when exceeding this core count, at least in the current implementation, concerns the mixing of input and output density matrices in the self-consistency loop of DFTB+, which is not yet MPI parallelized. This also affects the memory consumption. For some mixers (e.g. modified Broyden’s method [56]) the history of all previous $\Delta P_{\mu\nu}(\bm{g})$ is kept in storage, which becomes unfeasible for extremely dense k-point samplings or an unusual large number of self-consistent steps to reach convergence. Fortunately, mixing schemes with limited memory (e.g. modified Anderson’s method [57]) are readily available within DFTB+.

IV Convergence behavior

IV.1 Polyacene series

The important class of $\pi$-conjugated polymers has spawned numerous successful candidates for devices like organic light-emitting diodes (OLEDs) [58, 59, 60], organic field-effect transistors (OFETs) [61, 62, 63], polymer solar cells (PSCs) [64, 65, 66] and the growing field of organic electronics in general. One representative of this class is the \cfC_4n+2H_2n+4 series, forming the polyacene oligomers. We choose this linear molecular chain due to its relevance as a previous benchmark system for range-separated DFT. For many $\pi$-conjugated polymers (semi-)local DFT fails to describe the bond length alternation (BLA) and band-gap correctly. While Hartree-Fock (HF) overestimates BLA significantly, (semi-)local DFT is known to underestimate it. In fact, Körzdörfer et al. [67] suggested that the many-electron self-interaction error (MSIE) of HF and DFT approaches correlate with the BLA error and MSIE minimization is key (but not the only issue) to obtaining accurate BLAs. These quasi one-dimensional systems with low environmental screening also provide a stringent test for the removal of the divergence in the exchange interaction. In addition, the polyacenes feature a well-defined finite molecular limit, which can be used to verify the periodic RSH-DFTB implementation proposed here.

Since (semi-)local DFTB is directly derived from DFT, it also inherits its shortcomings. A particularly severe deficiency of PBE-DFT(B) is that the polyacene series becomes metallic for increasing chain length [68]. We employ the ob2-1-1 parameters [69], created for the purely long-range corrected LCY-BNL functional, to demonstrate that our implementation of periodic RSH-DFTB converges to the same, finite band-gap as the already available non-periodic formalism. We provide an estimate of the polyacene bandgap for the family of purely long-range corrected functionals, as calculated by the FHIaims code [70, 37, 38] on the LC-$\omega$PBE [71] level of theory using intermediate basis settings and a range-separation parameter of $\omega=0.3\,{a_{0}^{-1}}$, yielding $E_{\mathrm{gap}}=2.8\,\mathrm{eV}$ 111Due to technical difficulties experienced with the latest release version 221103 of the FHIaims code, this calculation has been carried out with version 210716.3 instead.. Since the specific LC-BNL functional is not yet available through FHIaims (release version 221103) 222FHIaims in version 221103 is interfaced with a recent version of the libXC library (5.1.7), that, in principle, provides an implementation of the BNL functional. However, LDA based hybrid functionals like BNL do not seem to be supported by FHIaims in version 221103. Furthermore, to the best of our knowledge, a manual adaption of the range-separation parameter $\omega$ is not yet supported, including the current development version of FHIaims, rendering a quantitative comparison of the ob2-1-1 parameters with the literature parametrization of the BNL functional [17], with $\omega=1\,a_{0}^{-1}$, meaningless., this value does not allow for a quantitative comparison with results obtained by the ob2-1-1 parameters, however, offers valuable guidance from first principles. Figure 3 illustrates the convergence of the polyacene band-gap convergence for the present periodic $\Gamma$-point and k-point implementation, in direct comparison to the non-periodic case. The calculations are based on the (unrelaxed) primitive unit cell of polyacene (k-point implementation), as listed in Structure S1 of the Supplemental Material [43], and supercells (the $\Gamma$-point implementation) built from it. In the case of the non-periodic implementation, the supercells are converted into clusters and properly passivated by additional hydrogen atoms at the chain ends (bond-length of C-H units: $1.1\,\text{\AA}$). The k-point sampling $1\!\times\!1\!\times\!(2\leq n\leq 100)$ is chosen according to the Monkhorst-Pack [74] scheme, where the polyacene chain is oriented along the $z$-direction and vacuum inserted in $x$- and $y$-direction.

What immediately stands out is that all implementations converge towards the same Kohn-Sham gap, which is an essential step in the validation 333Another strategy for verifying the correct implementation that has proven to be useful, is to test the translational invariance of the method by homogeneously shifting the coordinates (partly outside the unit cell) and verifying that the Kohn-Sham spectrum is unaltered. of the present method. A closer look at the convergence behavior of the general k-point implementations of the TCI and MIC schemes reveals strong fluctuations of the band-gap for non-convergent k-point samplings $1\!\times\!1\!\times n$, where $n\lesssim 20$. These fluctuations originate from artifacts in the band-structure and can be avoided by either manually reducing the Coulomb truncation cutoff in TCI or removing outer shells of unit cells inside the Wigner-Seitz cell of the BvK supercell in MIC. We take the opportunity of this specific case to highlight the pitfalls of periodic Fock exchange in RSH-DFT(B).

Figure 4 explicitly shows two fully converged band-structures of Figure 3, demonstrating that not only the gap size, but also all bands, calculated with TCI or MIC are virtually identical. Additionally, a non-convergent band-structure, originating from a density calculation with $1\!\times\!1\!\times 6$ k-point sampling is included. For this choice of parameters the band-gap collapses and individual bands exhibit an unphysical dispersion. Too small a BvK supercell does not allow for a natural decay of the density matrix, but rather introduces a spurious periodicity as described in Section III.1. The extent to which artifacts of non-convergent calculations based on coarse k-point samplings manifest themselves is system-specific. As an example shown in Figure S1 of the Supplemental Material [43], armchair graphene nanoribbons (AGNRs) turn out to be a much more benign system and convergence is achieved rapidly. We refer to Structures S3-S5 and Figures S2-S3 of the Supplemental Material [43] to obtain further investigations covering two-dimensional h-BN monolayer and GaAs bulk.

IV.2 Total energy convergence

The total energy is often considered to provide a solid indication of the convergence behavior of a system. In the literature [41, 39, 29] it is used for the sake of comparing implementations and to demonstrate convergence. However, this measure proved to be unreliable in many cases with regard to other properties, including band-gaps. We would therefore like to point out that, at least for periodic RSH-DFTB, further quantities of interest should also be considered when checking convergence.

Figure 5 illustrates that the absolute deviation in the total energy and band-gap with respect to their fully converged reference decreases at different rates, with slower convergence of the band-gap. This phenomenon appears to be independent of the dimensionality of the system and applies, e.g., for three-dimensional GaAs bulk and two-dimensional h-BN monolayer included in Figure 5. To obtain a fully converged total energy and band-structure as reference, we employed a GaAs density calculation with a $20\times\!20\times\!20$ and h-BN monolayer calculation with $20\times\!20\times\!1$ Monkhorst-Pack k-point sampling respectively.

V Polarization-induced gap renormalization

Renormalization of the fundamental band-gap in molecular crystals by electronic polarization [76] is of central importance for organic electronics [58, 59, 60, 61, 62, 63] and photovoltaics [64, 65, 66]. Going from the molecule in gas phase to a molecular crystal with relative dielectric constant $\varepsilon_{\infty}$ (orientationally averaged and ion-clamped) leads to shrinkage of the fundamental gap. This renders the resulting material well-suited for practical applications that require reduced optical gaps. Due to the electronic polarization of the crystalline dielectric medium, the energy required to create a quasi-hole is reduced compared to its molecular phase, whereas creating a quasi-electron releases more energy [77]. In other words, the ionization potential (IP) and electron affinity (EA) decrease and increase respectively.

Today’s standard repertoire of exchange-correlation functionals within DFT, including (semi-)local LDA/GGA as well as global and range-separated hybrids [15], do not properly treat long-range correlation effects and fail to describe the aforementioned gap renormalization, even qualitatively [78]. While many-body perturbation theory, especially Hedin’s GW approximation [79] to the electron’s self-energy, $\Sigma$, captures these renormalization effects, only recent screened range-separated hybrid functionals [77] include this effect at the considerably cheaper level of DFT.

According to the general CAM partitioning of the electron-electron interaction introduced by Eq. (1), the limiting behavior of the long-range part, that we treat in an exact Fock-like manner, is determined by the parameters $\alpha$ and $\beta$. For small distances as $r\to 0$, the $\alpha/r$ contribution prevails, whereas the limiting case of $r\to\infty$ scales like $(\alpha+\beta)/r$. In the gas phase, the correct $1/r$ asymptotic decay is obtained if the condition $\alpha+\beta=1$ is fulfilled. In fact this only represents the special case of $\varepsilon_{\infty}=1$ and a generalization to arbitrary dielectric environments with asymptotic potential of $1/(\varepsilon_{\infty}r)$ requires that $\alpha+\beta=1/\varepsilon_{\infty}$. For the gas-phase fundamental gap to coincide with the HOMO-LUMO gap of generalized Kohn-Sham RSH-DFT [80, 81], the Fock-like exchange term is required to be asymptotically correct [82, 83] and the range-separation parameter $\omega$ should be tuned to obey the ionization-potential (IP) theorem [21, 84, 85, 86]. We follow Refaely-Abramson et al. [77] by non-empirically determining $\omega$, through minimizing the function

[TABLE]

for the gas-phase. Here $\varepsilon_{\mathrm{H,n}}^{\omega,\alpha},\varepsilon_{\mathrm{H,a}}^{\omega,\alpha}$ denote the energies associated with the HOMO of the neutral (n) and anionic (a) systems, tuned to match the respective IP obtained from total energy differences as closely as possible. Additional figures outlining the optimization process are provided in Section S4 of the Supplemental Material [43]. For the sake of comparing our results with Ref. 77, we do not attempt to optimize $\alpha$ from first principles, but rather chose $\alpha=0.2$, which proved to yield satisfactory results for small organic molecules. The employed exchange-correlation functional is CAMY-PBEh [87] (with Y indicating the range-separation function is of Yukawa type). In order to consistently compare the results obtained with our RSH-DFTB method, we resort to the molecular and crystalline geometries of Ref. 77, covering the prototypical conjugated molecules benzene and pentacene. The same reasoning applies to the choice of the scalar dielectric constant, which is also taken from the cited reference. A full overview of the resulting functional parameterization, including the optimally tuned range-separation parameters $\omega$, is provided in Table S2 of the Supplemental Material [43]. To avoid time-consuming re-parameterization, the generated Slater-Koster files are based on the ob2-1-1 [69] parameters, which are expected to perform well in combination with long-range corrected functionals. Additionally, results obtained by a slight modification to the ob2-1-1 parameters are shown, this decreases the density compression radius of the carbon species to a value of $r^{\mathrm{dens}}_{\mathrm{C}}=3.3\,a_{0}$. We refer to Table S1 of the Supplemental Material [43] for a detailed listing of the electronic parameterization. Figure 6 compares the fundamental gaps obtained using different levels of theory.

The PBE-DFTB calculations were carried out based on the pbc-0-3 parameters [89]. Since the hydrogen and carbon electronic structure of pbc-0-3 is virtually identical to mio-1-1 [90]—the usual choice for biological or organic molecules within second-order DFTB—only one of the two options is included for comparison. In line with expectations, PBE-DFTB is unable to capture changes in the surrounding dielectric medium and severely underestimates not only the molecular but even the crystalline fundamental gap. In stark contrast, optimally tuned screened range-separated hybrid (OT-SRSH) DFTB as derived and implemented in this work captures the renormalization phenomenon at least qualitatively. However, the gas-phase HOMO is slightly too low in both systems, which also affects the bulk phase due to the chosen crystalline gap alignment. To maintain comparability with Ref. 77, the crystalline gap was aligned symmetrically around the middle of the respective gas-phase gap. Quantitative agreement would require a more in depth re-parameterization, which is outside the scope of this work. We emphasize that OT-SRSH-DFTB is systematically extending the domain of accessible exchange-correlation functionals within (periodic) DFTB by only introducing system-specific adjustable parameters. These are determined from first principles in a well-founded tuning process, rather than being subject to empirical fitting.

VI Summary and outlook

We have derived, implemented and tested Hartree-Fock exchange in the density functional tight-binding (DFTB) method for periodic systems beyond the $\Gamma$-point. By applying the usual DFTB approximations to matrix elements of the Fock exchange operator, we arrived at a real-space formulation of the exchange Hamiltonian that allows for efficient implementation in the DFTB+ code. To avoid artifacts owing to the artificial Born–von Kármán periodicity of the density matrix, we resort to either a truncated Coulomb interaction or a minimum image convention. These two methods proved to converge to the same limit, when coupled to accurate Brillouin zone samplings. Scaling with system size and parallel performance of the developed routines indicate that systems with thousands of atoms are well within reach of RSH-DFTB and that computational resources are exploited efficiently. Pre-screening of products of density matrix and overlaps, in combination with an iterative construction of the Hamiltonian allows us to further reduce the computational cost. Convergence behavior and the pitfalls of periodic RSH-DFTB are demonstrated for the polyacene series, showing that periodic and non-periodic implementations converge towards the same limits, and that the total energy is often not a reliable indicator for the convergence of other properties like the band-gap. In line with $G_{0}W_{0}$ calculations, screened RSH-DFTB shows the correct polarization-induced gap renormalization in benzene and pentacene molecular crystals, at a heavily reduced computational cost compared to first principles methods.

An in-depth benchmarking of the accuracy of periodic RSH-DFTB is currently subject of ongoing investigations and will be the topic of a future article.

Acknowledgements.

T. v. d. H.  and T. A. N. acknowledge financial support from the German Research Foundation (DFG) through Grant No. FR2833/76-1. The simulations were performed on the HPC cluster Aether at the University of Bremen, financed by the German Research Foundation (DFG) within the scope of Zukunftskonzept 66 “Ambitioniert und agil”, Bremen (GZ ZUK 66/1-2015).

Data availability

The data that supports the findings of this study are available within the article and its Supplemental Material [43].

Appendix A Hubbard parameter of RSH-DFTB

In analogy to conventional DFTB, the Hubbard $U$ of RSH-DFTB is derived by requiring its equality with RSH-DFT. Further, the xc-kernel kernel $f_{\mathrm{xc,loc}}$ is assumed to vanish [24] for off-site elements $A\neq B$, whereas the on-site elements $A=B$ cover the full exchange-correlation contributions and read as

[TABLE]

Note that, by definition, the function $\gamma^{\mathrm{fr,HF}}=\lim_{\omega\to\infty}\gamma^{\mathrm{lr,HF}}$ is always contained in the results for the more general screened Coulomb kernel. Since the Hubbard $U$ is obtained from a single atom RSH-DFT calculation, the following conditions apply

[TABLE]

By utilizing Eqs. (52) to (55), the total RSH-DFTB energy of the single atom in terms of occupation numbers $n_{i}$ is [24]

[TABLE]

with terms linear in $n_{i}$ indicated by $\mathcal{O}(n_{i})$. For the highest occupied shell it holds that for $n$ electrons, equally distributed over the shell, give an orbital occupation $n_{i}=n/d_{l}$, where the shell degeneracy is $d_{l}=2l+1$. Inserting $d_{l}$ into Eq. (57), yields

[TABLE]

Calculating the second derivative of Eq. (58) with respect to the shell occupation $n$ then becomes straightforward:

[TABLE]

In practice, the Hubbard $U$ of RSH-DFT is obtained by numerically calculating the second derivative of the eigenvalue $\varepsilon_{\mathrm{HOAO}}^{A}$ of the highest occupied atomic orbital of species $A$ with respect to its occupation $n_{\mathrm{HOAO}}$

[TABLE]

By enforcing $U^{\mathrm{RSH-DFTB}}\overset{!}{=}U^{\mathrm{RSH-DFT}}$, we finally end up with the Hubbard $U$ of RSH-DFTB

[TABLE]

an equation that should be solved numerically to obtain the decay constant $\tau_{A}$.

Appendix B $\mathrm{\Gamma}$-point approximation

For periodic RSH-DFTB calculations that are restricted to only the $\Gamma$-point, we resort to the TCI scheme exclusively. In the $\Gamma$-point-only approximation, any phase factors vanish and the $\bm{g}$-summation of Eq. (III.1) is carried out in advance, yielding new $\gamma$-integrals

[TABLE]

With $\gamma_{MN}^{\mathrm{TC}}$ truncated spherically and the $\bm{g}$-summation covering the entire crystal, the value of $\tilde{\gamma}_{MN}^{\mathrm{TC}}$ eventually becomes independent of the shift, $\bm{l}$, which greatly facilitates the pre-tabulation of $\tilde{\gamma}_{MN}^{\mathrm{TC}}$ for all element combinations.

From this we infer that the exchange Hamiltonian at the $\Gamma$-point $\Delta{\bf H}^{x,\mathrm{lr}}(\bm{\Gamma})$ can be constructed from matrix-matrix multiplications of dense overlap and density matrices

[TABLE]

with $\tilde{\gamma}^{\mathrm{TC}}$ denoting a super-matrix with the same shape as $S(\bm{\Gamma})$ and $\Delta P(\bm{\Gamma})$. Element-wise multiplication is indicated by $\odot$, i.e. the Hadamard product. Eq. (B) is a direct generalization of the non-periodic algorithm that is already implemented in the DFTB+ [25] code. Establishing $\tilde{\gamma}^{\mathrm{TC}}$ is straightforward as all elements of the diatomic block between atoms $M$ and $N$, containing orbitals {$\mu[M],\nu[N]$}, takes the same value: $\tilde{\gamma}_{MN}^{\mathrm{TC}}$, i.e. $\tilde{\gamma}^{\mathrm{TC}}_{\mu[M]\nu[N]}=\tilde{\gamma}_{MN}^{\mathrm{TC}}$. This algorithmic solution is appealing, since it is exact in the sense that no integral screening, as described in Section III.4, is required and since dense matrix-matrix multiplications can be performed efficiently and are easy to parallelize. A neighbor-list based algorithm for $\Gamma$-point calculations that includes integral pre-screening is implemented as well. However, for the systems so far tested, we observed a significantly higher efficiency for the matrix-multiplication based algorithm, therefore we refrain from discussing the list-based approach in detail.

Analogously to Eq. (B), the total energy and its gradients, i.e. atomic forces, can be expressed in a similar fashion

[TABLE]

with Cartesian coordinates $R^{\alpha}_{M}$ of atom $M$, $\alpha\in\{x,y,z\}$ and ${\bf A}^{\text{sym}}=\frac{1}{2}({\bf A}+{\bf A}^{\top})$ denoting the symmetric part of matrix ${\bf A}$.