Robust Periodic Fock Exchange with Atom-Centered Gaussian Basis

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Robust Periodic Fock Exchange with Atom-Centered Gaussian Basis Sets Andreas Irmler,† Asbjörn M. Burow,‡ and Fabian Pauly*,§,† †

Department of Physics, University of Konstanz, Universitätsstraße 10, D-78464 Konstanz, Germany Department of Chemistry, Ludwig-Maximilians-Universität (LMU) Munich, Butenandtstraße 7, D-81377 Munich, Germany § Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0395, Japan

Downloaded via KAOHSIUNG MEDICAL UNIV on August 25, 2018 at 06:19:50 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: In this work, we present a robust implementation of the periodic Fock exchange for atom-centered Gaussian-type orbitals (GTOs). We discuss the divergence, appearing in the formulation of the periodic Fock exchange in the case of a finite number of k-points, and compare two schemes that remove it. These are the minimum image convention (MIC) and the truncated Coulomb interaction (TCI) that we use here in combination with k-meshes. We observe artifacts in four-center integrals of GTOs, when evaluated in the TCI scheme. They carry over to the exchange and density matrices of Hartree−Fock calculations for TCI but are absent in MIC. At semiconducting and insulating systems, we show that both MIC and TCI yield the same energies for a sufficiently large supercell or k-mesh, but the self-consistent field algorithm is more stable for MIC. We therefore conclude that the MIC is superior to TCI and validate our implementation by comparing not only to other GTO-based calculations but also by demonstrating excellent agreement with results of plane-wave codes for sufficiently large Gaussian basis sets. space, with pioneering work by Gygi and Baldereschi,21 later by Spencer and Alavi,22 and finally by Sundararaman and Arias.23 Approaches for building the Fock exchange matrix in real space are the corresponding reformulation of the reciprocal space description of the truncated Coulomb interaction (TCI) by Guidon et al.19 as well as the minimum image convention (MIC) by Tymczak et al.20 Another possibility to perform robust calculations of periodic Fock exchange is the use of screened Coulomb kernels, as introduced by Heyd et al.18,24,25 in the context of range-separated exchange-correlation functionals. This approach can greatly accelerate the time needed for building the Fock matrix by spatially confining the nonlocal exchange interactions. It is therefore particularly promising for such periodic systems, where large simulation cells are inevitable. While different schemes for preventing the divergence of periodic Fock exchange were compared at realistic materials for plane-wave approaches,23 such a systematic study is presently missing for atom-centered basis sets. Furthermore, MIC and TCI were developed and applied at the Γ-point only for Gaussian-type orbitals (GTOs).19,20 Our work fills this gap by extending MIC and TCI implementations for a GTO-based code to arbitrary crystal structures and k-meshes. Finally, we compare both truncation schemes for a variety of crystals in order to find an optimal implementation of exchange in periodic systems.

1. INTRODUCTION Fock exchange is a crucial ingredient for accurate descriptions of the electronic structure in periodic systems. Despite known deficiencies, for instance, related to electronic self-interaction, the treatment of charged states and related band gap underestimation, or the proper characterization of spin states, density functional theory (DFT) presently remains the work horse of electronic structure calculations for solid-state systems.1,2 While pure DFT exchange correlation functionals hardly improved over the last years, so-called hybrid functionals became more and more popular, which are constructed by mixing Fock exchange into local and semilocal exchangecorrelation functionals.3−5 In addition, an active research area is the implementation of periodic boundary conditions for correlation methods such as the MP2 scheme of Møller−Plesset perturbation theory,6,7 the coupled cluster ansatz8 or the random phase,9,10 and GW11 approximations, which were presented recently. All of these methods use in some way the Fock exchange. Pioneering work on Fock exchange was performed by Pisani et al.,12 who made periodic Hartree−Fock (HF) calculations publicly available in the program package CRYSTAL early on. Nowadays, periodic HF implementations are available for planewave-based codes13−16 and also for software using localized numerical17 and Gaussian-type18−20 basis functions. A proper formulation of periodic Fock exchange is not straightforward because of the intrinsic divergence for any finite supercell size or equivalently any finite k-mesh. Different schemes exist for treating the integrable singularity in reciprocal © XXXX American Chemical Society

Received: February 5, 2018 Published: August 6, 2018 A

DOI: 10.1021/acs.jctc.8b00122 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation In this paper, we thus present the implementation of periodic Fock exchange using unscreened and screened Coulomb interactions. The routines developed are included in the periodic electronic structure module RIPER26 of the TURBOMOLE program suite,27 which uses GTOs. Our work is structured as follows. In Section 2, we derive the periodic HF expressions as we use them in our implementation. We focus on calculations with finite k-meshes leading to a finite supercell size. By comparing the two truncation schemes MIC and TCI in Section 3, we identify the optimal periodic exchange implementation regarding stability in the self-consistent field (SCF) method and convergence of total energies with supercell size. Finally, we discuss benchmark calculations in Section 4, where we compare our results with published data of other groups. In this context, we present, for instance, lattice constants, bulk moduli, and band transitions for selected insulating and semiconducting materials, which are in excellent agreement with previously published GTO and projector augmented wave (PAW) results. We close this work with performance results given in Section 5 and with a summary of our conclusions in Section 6.

∫Ω d3r Ψik*(r)Ψkj ′(r) = δijδ kk′ leading to independent normalization conditions

∑ Cμki*Sμνk Cνkj = δij for the orbital expansion coefficients overlap matrix in reciprocal space k Sμν =

l Sμν =

∑ e iklϕμl(r) l

EX = −

+ L) =

1 Nk

∑ Sμνk e−ikl (10)

k

1 Nk

∑ ∑ fik f jk′ ∫ ij

Ω

kk′

d3r1

∫ d3r2

Ψik *(r1)Ψkj ′(r1)Ψkj ′*(r2)Ψik (r2) |r1 − r2|

(11)

where i,j run over all bands and k,k′ run over the set defined in eq 4. The occupation numbers f ki range from zero to one. Spatial integration over r1 is restricted to the BvK supercell, while integration over position r2 is conducted over the entire real space. Inserting GTO-based orbital expansions, eq 11 can be rewritten as32

(1)

EX =

1 2

∑ Dμνl Kμνl (12)

μν l

in terms of the real space Fock exchange matrix l Kμν =−

(3)

1 2

∑ ∑ Dκλn − m(μ0 κ m|λ nν l) (13)

κλ mn

and the real space density matrix Dlμν obtained by Fourier backtransformation, cf. eq 10, of k Dμν =

∑ 2fik CμkiCνki*

(14)

i

We use here the Mulliken notation of the four-center electronrepulsion integrals (ERIs) (μ0 κ m|λ nν l) =

by introducing reciprocal primitive lattice vectors bj and Born− von Karman (BvK) boundary conditions31 Ψik (r)

(9)

For closed shells and crystal orbitals normalized to the BvK supercell, the Fock exchange energy of a single unit cell is defined as

employing real-valued GTOs ϕlμ(r) = ϕμ(r − Rμ − l) centered at atom positions Rμ in all unit cells l. Consequently, each sum over primitive lattice translation vectors shown in this subsection runs over the entire Bravais lattice. The wavevectors within the first Brillouin zone may be restricted to Nk = NUC different vectors30 mj k= ∑ bj , mj = 0, ..., Nj − 1 Nj j = 1,2,3 (4)

Ψik (r

(8)

∫ d3rϕμ0(r)ϕνl(r)

l Sμν =

(2)

1 NUC

∑ e iklSμνl

where superscript 0 stands for the reference unit cell and integration is conducted over the entire real space. The inverse transformation is performed with

expanded in terms of Bloch basis functions ϕμk (r) =

on every k-point. The

is obtained by a discrete Fourier transformation of the real-space overlap matrix

∑ Cμkiϕμk (r) μ

Ckμi

l

with integer numbers Nj. We use crystal orbitals belonging to bands i and wavevectors k, Ψik (r) =

(7)

μν

2. THEORY 2.1. Self-Consistent Field Equations and Fock Exchange. We start with a summary of the HF and Kohn−Sham (KS) equations for periodic systems, employing atom-centered GTOs. While our routines work both for open-shell and closedshell systems, we use in the following for simplicity the restricted closed-shell formulation.12,28,29 The periodic system or infinite ideal crystal is modeled by repetitions of unit cells l = ∑j=1,2,3ljaj with integers lj and primitive cell vectors aj. A supercell L=∑j=1,2,3ljLj consists of a finite number of NUC = N1N2N3 unit cells spanned by the supercell vectors Lj = Nja j , j = 1,2,3

(6)

∬ d3r1d3r2ϕμ0(r1)ϕκm(r1)g(r12)ϕλn(r2)ϕνl(r2) (15) −1

r−1 12 .

with the Coulomb metric g(r12) = |r1 − r2| = From minimization of the total energy, the canonical HF or KS crystal orbitals belonging to wavevector k are determined as solutions of the corresponding Roothaan−Hall equations

(5)

that constrain the Bloch orbitals to periodic functions on the superlattice defined by the base vectors Lj. The crystal orbitals are orthonormal within the volume Ω = |L1(L2 × L3)| of the BvK supercell

∑ Fμνk Cνki = ϵik ∑ Sμνk Cνki ν

B

ν

(16) DOI: 10.1021/acs.jctc.8b00122 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation with orbital energies ϵki and the Fock matrix in reciprocal space Fkμν obtained from the corresponding matrix in real space l l l l Fμν = Tμν + Jμν + Y μν

11 diverges due to the long-range Coulomb operator. This divergent behavior transfers to the long-range lattice sum within the real space expressions, as we demonstrate now. If the BvK supercell is finite, then the crystal orbitals and thus the density matrix are invariant for any superlattice translation vector L

(17)

by a transformation analogous to eq 8. The evaluation of the kinetic energy matrix Tlμν is identical to the molecular case, and those of the Coulomb matrix Jlμν follows ref 26. In eq 17, Ylμν represents a method-specific linear combination of exchangel correlation matrices Xμν from local or semilocal density functional approximations, which are calculated as published in ref 33, and of the nonlocal Fock exchange matrix Klμν, defined in eq 13. The evaluation of Klμν is the subject of the present work. Though explicitly shown here only for Fock exchange, see eq 13, the Coulomb and local or semilocal exchange-correlation contributions to the Fock matrix also depend on the density matrix. Equation 16 is thus solved iteratively until selfconsistency is reached. By setting Ylμν = (1−α)Xlμν + αKlμν in eq 17, the mixing parameter is α = 0 for local or semilocal DFT, α = 1 for HF, and 0 < α < 1 for hybrid density functionals. The range-separated hybrid functional approximation of Heyd, Scuseria, and Ernzerhof (HSE) is constructed in a similar fashion by replacing a fraction β of the short-range (SR) contribution to the semilocal exchange XPBE,SR (ω) by the same fraction of SR Fock exchange X KSR(ω)24 Y

HSE

(ω) = X

PBE

−

βXXPBE,SR (ω)

SR

+ βK (ω)

l l+L Dμν = Dμν

(19)

In eqs 12 and 13, we may therefore split up each lattice sum into a sum over all translation vectors of the Bravais superlattice and a sum over Nk primitive translation vectors within a single BvK supercell. If we indicate the latter by BvK above the summation symbol, we can write EX =

1 2

BvK

∑ ∑ Dμνl K̃μνl μν

(20)

l

with l K̃ μν =

1 2

∑ Kμνl + L = − ∑ DpqCμ0p,ν lq L

(21)

pq

Here, we have summarized indices in tuples p = (κ,m,M) and q = (λ,n,N), yielding brief notations for the density matrix n−m n + N − (m + M) Dpq = Dκλ = Dκλ

(18)

(22)

and for lattice sums of ERIs

where indices μ,ν,l are suppressed and the semilocal density functional approximation of Perdew, Burke, and Ernzerhof (PBE) is employed.34 The SR Fock exchange matrix KSR(ω) is directly obtained from eq 13 by using g(r12) = erfc(ωr12)/r12 in the underlying eq 15 with the attenuation parameter ω. In this paper, we employ two different parametrizations for the HSE functional. The mixing parameter is in both cases β = 0.25, but we use different values for ω. Krukau et al.35 proposed the use of the parameter ω = 0.11 a.u.−1 in eq 18. We refer to this parametrization as HSE06. In order to compare our results with the publication of Heyd and Scuseria,18 we also use the parameters ω = 0.15/ 2 a.u.−1 for the Fock exchange matrix KSR(ω) and ω = 21/3·0.15 a.u.−1 for the PBE part XPBE,SR(ω) (see ref 36). Results with this parametrization are named HSE03 below. 2.2. Finite Born−von Karman Supercells. A Bloch basis restricted to a finite BvK supercell is not complete with respect to all possible quantum numbers k of the infinite crystal and leads to finite-size errors. A complete Bloch basis can be approached by taking the limit Nj → ∞ for j = 1,2,3 so that sums over wavevectors turn into integrals, e.g., in eqs 10 and 11.31 For practical calculations of real solids, it is common to evaluate expressions with finite BvK supercells or equivalently on a finite grid of k-points. A reasonable choice of the supercell size for accurate results depends on the type of solid, the investigated property, and the level of electronic structure description.37 However, exchange interactions of finite BvK supercells contain undetermined contributions, and their evaluation thus requires additional corrections, as we explain in the following. As pointed out in previous work,21,38 the undetermined contributions can be most fundamentally understood by considering the Fock exchange energy in the orbital basis, eq 11, where orbital products Ψik*(r)Ψjk′(r) become latticeperiodic and real-valued functions of space with nonzero charge if and only if i = j and k = k′. In this case, the integral over r2 in eq

Cμ0p , ν lq =

∑ (ρμ00p |ρνLlq )

(23)

L

between real-valued charge distributions of GTO products ρμ00p = ϕμ0ϕκm + M

(24)

ρνLlq = ϕνl + Lϕλn + N + L

(25)

These GTO products contribute to the exchange matrix only for sets of tuples p and q within a short-range of the functions ϕ0μ and ϕl+L ν , respectively, where the cells m and n of p and q, respectively, are always within the BvK supercell. While the set of vectors l is also restricted to a single BvK supercell, the sum over the superlattice vector L in eq 23 is long-ranged and therefore converges only under well-known conditions.39 Divergent L-sums cannot be avoided, if both GTO products in eqs 24 and 25 are charged, unless eq 23 is modified in such a way that the exchange determined for the infinite BvK supercell is virtually recovered for finite supercells of moderate or large size. In the following Sections 2.2.1 and 2.2.2, we revisit two previously published modification strategies. The modifications can be motivated by first investigating the density matrix elements Dpq in the BvK supercell. It is known that the values of all Dpq drop to zero for large distances between p and q, where p is bound to the local region around the primitive reference cell.37,40−43 However, the BvK conditions in eq 5 lead to the unphysical periodic behavior of Dpq shown in eq 19. For infinitely large BvK supercells, corresponding to Nj → ∞ for j = 1,2,3, all density matrix elements drop to zero toward the BvK supercell boundary. In this case, the use of the BvK superlattice is an artifact, and all supercell translations may be neglected. In this way, eq 23 yields the well-defined value (ρ0μ0p|ρ0νlq), and the lattice sum over the primitive vectors l in eq 20 becomes longranged. But at the same time, the sum is damped out by the density matrix, yielding a corresponding well-defined value for C


Article

Journal of Chemical Theory and Computation the Fock exchange energy.23 For finite BvK supercells, the general idea is thus to suppress the interactions between the BvK images. As a side note let us remark that an equivalent formulation that avoids reference to a particular basis set can be obtained b y c o n s i d e r i n g th e o n e - e l e c t r o n d e n s i t y m a t r i x ρ(r1, r2) = ∑i,k2f ki Ψki (r1)Ψki *(r2). From this definition and eq 5, the periodicity ρ(r1,r2 + L) = ρ(r1,r2) follows for a finite number of employed k-points, while ρ(r1,r2) should decay with distance |r1−r2| in the real crystal.37,40−43 The exchange energy per unit cell in eq 11 can now be expressed as EX = −1/(4Nk)∫ Ωd3r1∫ d3r2ρ(r1,r2)ρ(r2,r1)/|r1 − r2|. While the unphysical periodicity of the offdiagonal elements of the oneelectron density matrix ρ(r1,r2) will impact the exchange energy of a crystal, described by a finite BvK supercell, this problem does not appear for the Coulomb energy, EJ = 1/(2Nk)∫ Ωd3r1∫ d3r2ρ(r1,r1)ρ(r2,r2)/|r1−r2|. EJ depends only on the diagonal part ρ(r1,r1), i.e., the electron density, which is naturally periodic on every single unit cell. 2.2.1. Truncated Coulomb Interaction. Replacement of the long-range Coulomb metric in eq 23 by the TCI l o1/r12 if r12 < R c g TCI(r12) = o m o o else n0

the values of the function G0(ρ , T ) = 2πF0(T )/ρ

(34)

with the Boys function F0(T) are used to evaluate four-center ERIs of the Coulomb metric g(r12).44 For the gTCI(r12) metric, this function becomes19 2π F0(T ) ρ

G0(ρ , T , R c′) =

+

π 3/2 erf(R c′ − 2ρ

T ) − erf(R c′ + T

T) (35)

where R c′ = ρ R c is the scaled truncation radius. Higher order functions Gn(ρ,T,Rc′), and similarly functions Fn(T), are obtained by Gn(ρ , T , R c′) = ( −∂T)n G0(ρ , T , R c′)

(36)

where the required maximum n depends on the atomic angular momentum quantum numbers at the four centers of the considered ERI. We find a compact algebraic representation of eq 36 by assuming a minimal Rc′ value such that erf(R c′ + T ) is virtually one for any T. Using this assumption, the last two terms in eq 35 can be summarized to

(26) 19,22

avoids divergent lattice sums from the beginning. In the expression, Rc is an adjustable cutoff radius. In all our calculations, this value is chosen such that interactions with neighboring BvK supercells are eliminated; i.e., Rc is the maximum radius of a sphere that fits into the BvK supercell. Efficient and simple evaluation of ERIs with the TCI metric is possible based on a conventional recurrence scheme,44 where the underlying modified Boys functions of higher order may be obtained from a sophisticated numerical evaluation scheme designed for general metrics and using precalculated data maps.19 We demonstrate now that a simple and robust algebraic scheme can be used to evaluate the modified Boys functions for the TCI ERIs, provided that the selected Rc value is not too small. This restriction is motivated by our knowledge that BvK supercell dimensions should not be chosen too small in order to obtain reasonable results for the exchange interactions, and we show in the present work that this assumption is indeed valid for practical calculations. We start from four-center ERIs, eq 15, of primitive Cartesian GTOs centered around A

2π [F0(T ) − A 0(T , R c′)] ρ

G̅0(ρ , T , R c′) =

(37)

The first term is the Boys function from eq 34, used for the Coulomb integrals. The second term modifies the Boys function A 0(T , R c′) =

π erfc(R c′ − 4 T

T) (38)

Each higher order function G̅ n(ρ,T,Rc′) is evaluated with the higher order Boys function Fn(T), which may be obtained by a recurrence formula,45 and the higher order function An(T,R′c), defined in analogy to eq 36. Evaluation of the derivative (−∂T)nA0(T,Rc′) leads to the algebraic expression Ä 1 ÅÅÅÅ ′ A n (T , R c ) = Å(2n − 1)A n − 1(T , R c′) 2T ÅÅÇ ÉÑ 2 Ñ 1 + e−(R c′− T ) Pn(T , R c′)ÑÑÑÑ ÑÖ (39) 2

2

ϕa(r) = (x − Ax )ax (y − A y )ay (z − Az )az e−α | r − A|

with bivariate polynomials

(27)

with basis function exponent α and atomic angular momentum quantum numbers ax, ay, and az. For the four-center ERIs, common centers between the GTO products ϕaϕb and ϕcϕd P = (α A + β B)/ξ

(28)

Q = (γ C + δ D)/η

(29)

P1(T , R c′) = −1

(40)

i R′ y Pn(T , R c′) = −jjjj c − 1zzzzPn − 1(T , R c′) − ∂TPn − 1(T , R c′) k T {

(41)

We have not found any recurrence relation for these polynomials; however, compact expressions for low orders can be explicitly derived in terms of the variables u = 1/(2T) and v = R c′/ T , e.g.,

are defined employing summarized exponents ξ=α+β

(30)

η=γ+δ

(31)

Together with additional definitions

P2(u , v) = v − 1

(42)

ρ = ξη/(ξ + η)

(32)

P3(u , v) = −(v − 1)2 + uv

(43)

T = ρ|P − Q|2

(33)

P4(u , v) = (v − 1)3 − 3(u + v − 1)uv

(44)

D


Article

Journal of Chemical Theory and Computation This finally means that the functions G̅ n(ρ,T,R′c) can be efficiently calculated order by order starting from eq 37. In comparison to ordinary Boys functions, additional computationally demanding terms, namely, the complementary error function in eq 38, the exponential function in eq 39, and the square root of T, are computed only once for all orders. The other required operations are additions and multiplications, and most of them are summarized in the polynomials Pn(T,R′c). In our code, ERIs of the TCI metric are currently implemented up to n = 14. The minimal Rc′ value is set to 6 because for this argument the error function erf(R c′ + T ) deviates only by less than 2.2 × 10−17 from 1. In addition, instabilities in the evaluation of higher orders of the modified Boys functions are circumvented by use of this minimal Rc′ value. 2.2.2. Minimum Image Convention. In order to avoid the divergent lattice sum in eq 23, Tymczak et al. suggested that the MIC can be applied to GTO-based periodic exchange interactions.20 Their scheme is formulated in the Γ-point approximation only, and we have adapted it here to the use with k-points and arbitrary crystal lattices, as explained in the following. The MIC scheme manipulates the sum over superlattice vectors L in eq 23. Only a single superlattice vector LMIC is selected, which minimizes the distance between the charge L distributions ρ0μ0p and ρνlq over all L. Equation 23 then turns into Cμ0p , ν lq = (ρμ00p |ρνLlqMIC )

since it will be retained for a BvK supercell of increasing size. In L contrast, the ERIs (ρ0μ0p|ρνlq ) with L ≠ LMIC are mainly contracted with unphysical periodic images of the density matrix and are therefore not present in eq 45. For sufficiently large BvK supercells, the MIC procedure ensures a symmetric use of the two density matrices Dlμν and Dn−m κλ in eq 20 since n−m will be in the same range as l. In summary, artificial exchange interactions between periodic BvK images are virtually excluded in eq 45 for large enough supercells. By selecting a single superlattice vector LMIC in eq 23, the exchange matrix in eq 21 as well as the exchange energy in eq 20 become well defined. In our implementation, we assign center positions P and Q + L to the charge distributions ρ0μ0p = ϕ0μϕm+M and ρLνlq = κ l+L n+N+L ϕν ϕλ in eqs 24 and 25. In the case that all GTOs are primitive, P and Q + L are defined through the Gaussian product theorem, similar to what is shown in eqs 28 and 29. If the charge distributions contain contracted GTOs, we select the most diffuse basis function as the relevant primitive of each contracted GTO and determine the center position of the whole distribution again from primitives. Consequently, Cμ0p,νlq in eq 45 is obtained from the contribution of the superlattice vector LMIC, which yields the minimum distance between P and Q + L. We note that the expression for the exchange matrix of ref 20 is recovered, if we omit all lowercase lattice vectors in eqs 21 and 45, representing the unit cells in the BvK supercells. However, our technical realization to determine LMIC differs from the scheme proposed in ref 20. Our algorithm computes minimum distances in Cartesian coordinates. This means that we work exclusively with the vectors P and Q + L. We search for LMIC at L = 0 and in the shell of surrounding BvK supercells. For large enough BvK supercells, LMIC will correspond to one of the these supercells. In contrast, the scheme of Tymczak et al.20 minimizes the connection vector P − (Q + L) in the corresponding fractional coordinates. The Cartesian coordinate P is connected to the fractional coordinate p by

(45)

The MIC procedure can be motivated as follows. The artificial periodicity property of the density matrix in eq 19 results in the divergent lattice sum over L in eq 23. Given the charge distribution ρ0μ0p, the Coulomb integrals with distributions ρLνlq for all of the different L contribute to Cμ0p,νlq in eq 23, as LMIC visualized in Figure 1. From these distributions, ρνlq is located

P=M×p

(46)

In the expression, the transformation M = (L1,L2,L3) consists of a matrix, formed by the supercell vectors Lj. If the crystal is described in terms of primitive cell vectors aj, which form a nonorthogonal basis set, M will not be unitary. In this case, a minimization of P − (Q + L) in Cartesian coordinates is not equivalent to a minimization in terms of fractional coordinates. For this reason, we identify LMIC in terms of Cartesian coordinates.

3. COMPARING MIC AND TCI APPROACHES We note that TCI and MIC are fundamentally different. In TCI, the Coulomb interaction is manipulated. In contrast, in MIC the full Coulomb interaction is taken into account, but a divergent sum over a superlattice vector is restricted to a single term. Both schemes are constructed such that artificial interactions between periodic images of the density matrix are basically eliminated, and the correct expression for the exchange energy is recovered in the limit of an infinitely extended BvK supercell. The following results are obtained with our implementation of MIC and TCI approaches explained in the previous section. First, we demonstrate that interactions between GTO products may exhibit artifacts when using the TCI metric if the distance between the product centers is close to the value of the truncation radius Rc. As an example, we start by comparing results for the ERI (s0s0|sdsd) of four primitive s-type GTOs, which is once evaluated using the Coulomb metric and once

Figure 1. Two-dimensional scheme, visualizing the charge distributions appearing in eq 23. ρ0μ0p is depicted as a blue s-type Gaussian. It is located in the central BvK supercell at position P, as indicated by the blue dot. The green s-type curves represent the distributions ρLνlq, located at Q + L, as shown by green dots. These distributions are repeated in every BvK supercell due to the summation over L in eq 23. The red line represents the Wigner−Seitz supercell, centered around P.

inside the “Wigner−Seitz supercell” centered around ρ0μ0p and L yields the term (ρ0μ0p|ρνlqMIC) of eq 45, which is contracted with density matrix elements belonging essentially to the “WignerSeitz supercell”. By “Wigner−Seitz supercell”, we mean the supercell obtained from the Wigner−Seitz construction using the supercell vectors Lj (see eq 1) instead of the conventional primitive cell vectors aj. Equation 45 identifies the relevant term E


Article

Journal of Chemical Theory and Computation

fine k-meshes, where the MIC approach converges slightly faster. For the TCI scheme, results for smaller k-meshes are not obtained for two reasons. For the two smallest k-meshes used for h-BN and LiH, the scaled cutoff radius, Rc′, was below the required minimum of reliable ERI evaluation within the approximation discussed in Section 2.2.1 due to a small corresponding BvK supercell. For the other missing results, the TCI approach did not lead to SCF convergence. In order to find an explanation for the SCF instability of the TCI method, we analyze the Fock exchange and density matrices of the h-BN sheet. For this purpose, we plot in Figure 3 the absolute values of Dlμν and K̃ lμν (see eq 20) for the MIC and TCI approaches as functions of the primitive lattice vector l in the 31 × 31 BvK supercell. We select the diffuse s-type GTO with exponent 0.268 a.u.−2 of the nitrogen atom from the pobTZVP basis set46 for indices μ and ν. In Figure 3, we compare the matrices between MIC and TCI for one HF-SCF cycle, starting from the same initial density matrix, obtained from a converged PBE calculation. We choose this strategy because it reveals the direct impact of the two different schemes on the next SCF cycle. The initial density matrix, shown in Figure 3(e), is used to construct the exchange matrices in Figure 3(a) and (c). These exchange matrices are included in the Fock matrices used to obtain new orbital coefficients by solving the Roothaan−Hall eq 16. The density matrices of the next SCF cycle, shown in Figure 3(b) and (d), are constructed subsequently from the orbital coefficients via eq 14. In Figure 3, the plots of absolute values of exchange matrix elements exhibit nearly spherical patterns, where the absolute values decay monotonically and smoothly with increasing radius from the primitive reference cell for the MIC approach. For the TCI method, artifacts appear near Rc, which are similar to those of the single ERI shown in Figure 2. Beyond this radius, TCI exchange matrix elements drop off abruptly. Note that an element of the form K̃ lss contains many different four-center ERIs contracted with density matrix elements and thus accumulates the effects of the TCI from all these ERIs. For primitive cells inside the sphere of radius Rc, which are not too close to Rc, the exchange matrix values from MIC and TCI agree with each other for the large BvK supercell used in Figure 3. This explains identical HF energies for 31 × 31 k-points from converged SCF calculations with both approaches, as shown in Table 1. In contrast to the exchange matrices in this example, the density matrices exhibit a hexagonal pattern, reminiscent of the crystal symmetry. In addition, density matrix elements of the next SCF cycle increase at the Rc boundary for the TCI procedure but not for the MIC method. This means that the effects of the TCI metric on the exchange matrix transfer into the density matrix of the next cycle.

using the TCI metric. The centers of the GTOs, 0 and d, yield a distance d between the GTO products. Figure 2 shows that the

Figure 2. Values of ERIs for the Coulomb metric (Rc = ∞) and the TCI metric (Rc = 35 a.u.) plotted against the distance d between the centers of bra and ket distributions. The inset shows a magnified plot for d ≈ Rc. Primitive GTOs with exponents of 0.1 a.u.−2 are used.

ERI exhibits a monotonic decay with d for both metrics. For the Coulomb metric, the interaction between the GTO products decreases as 1/d, while the TCI metric leads to an interaction localized within a domain bound to the Rc radius. Next, we replace one s-type GTO in this ERI by a p-type GTO yielding (p0x s0|sdsd), where d is now oriented along the x-axis. Except for a sign switch in the short-range, the Coulomb metric yields a monotonic 1/d2 decay of the ERI in the long-range. The absolute value of the ERI in the TCI metric decreases also rapidly with d; however, for d ≈ Rc, it increases to a local maximum value. This maximum appears at d = Rc since one lobe of the p-type GTO is outside of the interaction radius of the TCI metric, while the other lobe of opposite sign is inside. Therefore, the contributions of the two lobes do not largely cancel each other as it is the case for the Coulomb metric, for which the interaction radius is infinite. Similar artifacts may appear for ERIs of other types of GTOs and for other GTO positions if the distance between the centers of the bra and ket distributions is close to Rc. In the next step, we investigate, if these artifacts have consequences for practical HF calculations. Self-consistent HF energies of an atomically thin hexagonal boron nitride (h-BN) sheet and of lithium hydride (LiH) in the rocksalt structure are shown in Table 1. For both systems, the energies of the MIC and TCI schemes converge toward the same limit with increasingly

Table 1. Self-Consistent HF Total Energies per Primitive Cell in a.u. for Two-Dimensional h-BN with a Lattice Constant of 2.503 Å and for LiH in the Rocksalt Structure with a Lattice Constant of 4.084 Å19,a h-BN

LiH

k-mesh

MIC

TCI

k-mesh

MIC

TCI

5×5 7×7 13 × 13 17 × 17 21 × 21 31 × 31

−79.29057349 −79.29071336 −79.29071980 −79.29071986 −79.29071986 −79.29071987

− −b −c −c −79.29071990 −79.29071987

5×5×5 7×7×7 9×9×9 11 × 11 × 11 13 × 13 × 13 19 × 19 × 19

−8.06060281 −8.06058890 −8.06058834 −8.06058830 −8.06058829 −8.06058829

−b −b −8.06058764 −8.06058831 −8.06058829 −8.06058829

b

a

Primitive unit cells with two atoms are calculated for different k-meshes using the pob-TZVP basis set.46 bBvK supercell too small. cSCF unstable. F


Article


Figure 3. Fock exchange and density matrix elements |K̃ lμν| and |Dlμν| in a.u. after the first iteration for MIC and TCI approaches, plotted against the primitive cell vector l within the 31 × 31 BvK supercell of h-BN. Indices μ and ν belong to the diffuse s-type GTO of the N atom in the pob-TZVP basis set. (a) |K̃ lμν| and (b) |Dlμν| for MIC. (c) |K̃ lμν| and (d) |Dlμν| for TCI. (e) Initial density matrix elements |Dlμν| from a converged DFT calculation with the PBE exchange-correlation functional. (f) Profiles along the dashed lines indicated in panels (a)−(d).

We show here the spatial behavior of |K̃ lμν| and |Dlμν| only for the specific case that μ and ν are s-type GTOs, but similar results can be found for other kinds of GTOs. While the artificial truncation effects of the TCI on GTO-based ERIs are still seen when used together with quite large BvK supercells, these effects are much more pronounced on an absolute scale in smaller BvK supercells. As this yields an increase of the absolute values of elements in the exchange and density matrices at the Rc boundary, we hypothesize that these erroneous elements amplify each other to some degree, leading to the SCF instabilities that we observe for the h-BN sheet with 9 × 9 to 17 × 17 k-points. In Figure 4, we compare self-consistent Fock exchange and density matrices of the h-BN sheet between the MIC and TCI approaches using the 31 × 31 and 51 × 51 BvK supercells and the pob-TZVP basis set. The MIC approach leads again to a smooth decay of both |K̃ lμν| and |Dlμν| with increasing radius from the primitive reference cell, where the density matrix exhibits a much more widely extended hexagonal star-like pattern than the exchange matrix. The larger BvK supercell allows for a progressive decay of the matrix elements toward zero close to the BvK supercell boundary. Absolute values of matrix elements at the boundary reach a magnitude of 10−12 a.u. for the exchange and of 10−9 a.u. for the density in the 51 × 51 supercell, while they just decrease to approximately 10−9 and 10−7 a.u., respectively, in the 31 × 31 supercell. In the case of the TCI method, the same patterns occur in the center of the BvK supercells. However, boundary effects at the truncation radius Rc are visible even for the 51 × 51 supercell, though much less pronounced compared to the results for the 31 × 31 supercell.

While the TCI exchange matrix elements beyond Rc drop abruptly to zero, the magnitude of density matrix elements in the same region does not decrease much below 10−7 and 10−9 a.u. for the 31 × 31 and 51 × 51 BvK supercells, respectively. Density matrix elements beyond Rc thus remain of a similar absolute size in MIC and in TCI, whereas the magnitudes of exchange matrix elements differ. Sundararaman and Arias23 showed an exponential convergence of EX with respect to the number of k-points used by performing single-shot calculations with wave functions obtained from DFT. With the present implementation, we can go beyond their work by studying the convergence of selfconsistent total energies, as shown in Figure 5(a), where we use the MIC truncation scheme and plot the finite-size error E(N1/3 k ) − E(∞) of the total energy E as a function of the total 1/3 1/3 number of k-points Nk = N1/3 k × Nk × Nk . The data for the HF and PBE methods exhibit an exponential decay of the finite-size error with respect to the number of k-points. For HF, the decay has a larger exponent compared to PBE. This result can be linked to the stronger localization of the HF density matrices compared to those of DFT, as it is apparent from the plots in Figure 5(b). For the missing data points of the silicon HF calculations with the def-SVP basis set on the 5 × 5 × 5 and 7 × 7 × 7 k-meshes in Figure 5(a), we find SCF instabilities similar to our previous observations for h-BN in Table 1 calculated with the TCI scheme. Our HF tests on Si, employing TCI and the def-SVP basis set, yield instabilities for k-mesh sizes up to 11 × 11 × 11. The MIC instead already converges for k-meshes starting from 9 × 9 × 9, as can be inferred from Figure 5. We want to emphasize that we use very accurate thresholds for the total G


Article


Figure 4. Self-consistent Fock exchange (left) and density (right) matrix elements |K̃ lμν| and |Dlμν| in a.u. for MIC and TCI approaches, plotted against the primitive cell vector l within the 31 × 31 and 51 × 51 BvK supercells of h-BN. Indices μ and ν belong to the diffuse s-type GTO of the N atom in the pob-TZVP basis set.

Figure 5. (a) Self-consistent total energy differences |E(N1/3 k ) − E(∞)| per primitive cell for PBE and HF methods and (b) maximum absolute 0 l elements of the self-consistent density matrices |Dlμν| for N1/3 k = 25 ordered by atom distances between shells μ and ν in the BvK supercell. LiH is 19 calculated in the rocksalt structure with a lattice constant of 4.084 Å and is described with the basis set of ref 47. For Si, we use the diamond structure with a lattice constant of 5.430 Å13 and the pob-TZVP46 as well as def-SVP48 basis sets. E(∞) is approximated by the energy obtained with a 31 × 31 × 31 k-mesh since |E(25) − E(31)| < 4 × 10−10 a.u. for LiH and Si.

energy, requesting energy differences in two subsequent SCF steps to be below 10−9 a.u. Less tight thresholds result in

converged SCF calculations for smaller BvK supercells. Our observations of SCF instabilities might also be influenced by the H


Article


discuss the timings and number of evaluated integrals, which are specified in Table 2, in Section 5. In analogy to previous work, we have extracted lattice constants a0 and bulk moduli B0 by fitting the Murnaghan equation of state52 to total energy versus volume curves. Let us compare our results for these two quantities in LiH with those of HF calculations performed with the VASP code using PAWs.53 Their lattice constant of 4.111 Å fits very well to our value of 4.104 Å, and also the bulk modulus of 32.9 GPa from the RIPER calculation is close to the published value of 32 GPa of the VASP group.53 On the basis of our procedures to compute the Fock exchange, we have implemented different hybrid functionals, such as PBE0,4 and different range-separated functionals, such as HSE0324,36 and HSE06, 35 into the RIPER module of TURBOMOLE (see eqs 17 and 18). To validate the implementation, we compare results for the three materials silicon (Si), diamond (C), and magnesium oxide (MgO) to those reported in the study of Heyd and Scuseria18 with the program package GAUSSIAN. Therefore, we use the PBE and HSE03 functionals in combination with exactly the same basis sets employed in their work, namely, 6-31G*, for silicon as well as diamond and the basis set constructed by Catti et al.55 for MgO. The latter basis is in the following simply labeled “Catti”. Also the crystal structures are chosen consistently: We use the diamond structure for silicon as well as diamond and the rocksalt structure for MgO. Table 3 shows the comparison of the evaluated equilibrium lattice constants and bulk moduli. Deviations of the lattice constants are below 0.12%. Bulk moduli differ by no more than 1.8%. Beyond these mechanical properties, bandgaps for bulk Si, diamond, and MgO were reported by Heyd and Scuseria18 with PBE and HSE03 functionals. For comparison, we have recalculated the gaps with our periodic implementation using the optimized lattice constants from ref 18, which are listed in Table 3 in the column called “GAUSSIAN”. All gap results are in excellent agreement, and deviations are below 2.3% for all materials and functionals: For Si, we obtain 0.70 and 1.28 eV for PBE and HSE03, respectively, which are 0.01 and 0.03 eV below the previously published gaps. For diamond, our PBE bandgap of 4.07 eV is 0.03 eV below their value, while for the HSE03 functional our simulation yields a bandgap of 5.34 eV, which is 0.09 eV lower than their result. For MgO, the same bandgap value of 4.91 eV is calculated with both codes for the PBE functional, while our bandgap of 7.04 eV for the HSE03 functional is 0.05 eV below their value. Since the PAW approach is the de facto standard for calculations of periodic systems, comparison of our results with the PAW study of Paier et al.13,56 based on VASP is of special interest. The deviations of our structural properties in Table 3 increase from GAUSSIAN to VASP. The average deviations between RIPER and VASP for lattice constants and bulk moduli are however still less than 0.45% and 8.8%, respectively. To be precise, the deviations between RIPER and VASP are below 0.45% and 0.80% for a0 and B0 in silicon, 0.28% and 1.2% in diamond, and 0.36% and 8.8% with the Catti basis set in MgO, respectively, with all three functionals tested. This means that deviations in B0 for MgO are particularly large. For band transitions, a similar picture of large deviations for MgO is obtained between RIPER and VASP: For all three functionals used, deviations are below 9.9% for silicon, 0.43% for diamond, and 17.1% for MgO. Note that on an absolute scale band

specific choice of the SCF convergence acceleration scheme, which is in our case a Γ-point-only direct inversion in the iterative subspace algorithm,49 as described in ref 26. In summary, both MIC and TCI become unstable for small BvK supercells but converge exponentially with N1/3 k to the same energy if the supercells become large enough. In particular, for the h-BN and silicon crystals, we have demonstrated that the MIC scheme is more stable than TCI. We therefore use the MIC approach for all further Fock exchange calculations in the present work.

4. BENCHMARK CALCULATIONS Let us continue our analysis of LiH with the basis set of ref 47, which was specifically constructed along the lines of polarization-consistent basis sets. In Table 2, we list HF total energies Table 2. Self-Consistent HF Total Energies per Primitive Unit Cell in a.u. for LiH with Different k-Meshes and FourCenter Integral Neglect Thresholdsa Threshold (a.u.)

Energy (a.u.)

CPU time (s)

no. integrals (106)

5×5×5

10−8 10−10 10−12

−8.0645683470 −8.0645680566 −8.0645680575

284 1062 2069

314 868 1838

7×7×7

10−8 10−10 10−12

−8.0645521469 −8.0645518712 −8.0645518718

568 1946 4078

614 1878 4210

9×9×9

10−8 10−10 10−12

−8.0645517341 −8.0645514583 −8.0645514589

1142 4091 9315

933 3162 7511

11 × 11 × 11

10−8 10−10 10−12

−8.0645517224 −8.0645514466 −8.0645514472

1653 6534 17517

1223 4581 11537

k-mesh

a

We compute primitive cells with two atoms and employ the same basis set and crystal structure as Paier et al.47 in their GTO-based HF studies, i.e., the rocksalt structure with a lattice constant of 4.084 Å. In addition, timings for a single construction of the Fock exchange matrix are given in seconds, which should be compared with the number of non-negligible integral batches calculated.

per unit cell using different sizes of the BvK supercell and different four-center integral neglect thresholds. We find that both parameters influence the value of the total energy. In our calculations, we use the Schwarz inequality50 as integral neglect criterion for the four-center ERIs. We emphasize that the Schwarz cutoff is the only adjustable threshold in our construction of the Fock exchange matrix. On the basis of the results for the different k-meshes and different four-center integral neglect thresholds, we can identify the converged HF energy per unit cell for the given basis set with a precision of better than 10−6 a.u. Our energy value is in good agreement with those published in the joint work47 of developers of CP2K and GAUSSIAN of −8.064545 a.u. The CRYSTAL group was also able to reproduce this result with their implementation.51 In order to avoid any effect of the integral screening, we employ a four-center integral neglect threshold of 10−12 a.u. for all of the following results. Furthermore, we combine primitive cells with a 25 × 25 × 25 k-mesh if not noted otherwise. We I


Article

Journal of Chemical Theory and Computation Table 3. Lattice Constants a0 and Bulk Moduli B0 for Different Materials and Functionalsa RIPER PBE

PBE0

Si a0 (Å) B0 (GPa)

5.491 88.5

5.455 99.4

C a0 (Å) B0 (GPa)

3.584 426

3.557 464

MgO a0 (Å) B0 (GPa)

4.243 162

4.197 177

MgO a0 (Å) B0 (GPa)

4.267 147

4.216 164

GAUSSIAN HSE06

HSE03

PBE

6-31G*

VASP

HSE03

PBE

PBE0

HSE06

5.454 99.1

5.469 87.8

PAW 5.433 99.0

5.435 97.7

3.557 458

3.574 431

PAW 3.549 467

3.549 467

4.200 175

4.258 149

PAW 4.211 169

4.210 169

− −

4.258 149

PAW 4.211 169

4.210 169

6-31G* 5.459 98.3

5.452 99.3

5.490 88.6

3.558 463

3.556 466

3.583 422

4.198 176

4.195 177

4.242 161

4.217 164

− −

− −

6-31G*

6-31G*

Catti

Catti

Catti-mod

Catti-mod

a Our RIPER results are calculated with the 6-31G* basis set54 for silicon as well as diamond. For MgO, they are obtained with the basis set presented by Catti et al.,55 called “Catti”, and a modified one, “Catti-mod”, described in the text. The calculations are compared with published values from Heyd and Scuseria,18 obtained with GAUSSIAN, who used the identical basis sets except for the newly added “Catti-mod”. Furthermore, we list PAW results published by Paier et al.13,56 and computed with VASP.

approaches are incomplete GTO basis sets. Since we can virtually reproduce the GAUSSIAN results, we examine the assumption of Paier et al.13 by calculations of bulk MgO employing an improved GTO basis set. Figure 6(a) shows a comparison between band structures, calculated for the PBE functional in RIPER with the Catti basis set and in QUANTUM ESPRESSO14 with PAWs.57 We observe larger deviations, as anticipated from the results shown in Table 4, which originate mainly from the unoccupied conduction bands. The discrepancies in the band structure can be decreased drastically by slightly modifying the basis functions on the Mg atoms, while keeping them unchanged on the O atoms. By adding a p-type and an s-type basis function on the Mg atoms, both with an exponent of 0.12 a.u.−2, and changing the exponent of the d-type basis function from 0.65 to 0.25 a.u.−2, we arrive at the modified Catti basis set “Catti-mod”. As a consequence of the adjustments, the bandstructures of RIPER and QUANTUM ESPRESSO now show a good overall agreement, as can be seen in Figure 6(b). Consequently, we find that deviations for the given band transitions in Table 4 typically shrink by almost an order of magnitude and are less than 1.1% for the transitions Γv → Γc and Γv → Lc. For the Γv → Xc transition, deviations reduce only to below 3.2%. This can also be seen in Figure 6(b), where the differences between bandstructures of RIPER and QUANTUM ESPRESSO are largest for the lowest conduction band around the X-point. Our results thus confirm that the deviations between PAW and GTO results indeed arise from a poor description of the electronic structure of the Mg atom in the Catti basis set, which can be substantially improved upon straightforward basis set modification. Let us finally point out that the adjustments in the basis set Catti-mod do not only influence the conduction bands in Figure 6 but also modify lattice constants and bulk moduli in Table 3, which we compare for all three functionals between RIPER and VASP. While values of a0 for the Catti basis set were smaller than those of VASP, they are now bigger for Catti-mod. The opposite trend of a larger value for Catti and a smaller one for Catti-mod holds for B0. Overall, the agreement between VASP and RIPER is slightly improved: While the deviations for the Catti results

transitions for silicon and diamond agree by better than 70 meV, while they deviate by up to 1.54 eV for MgO. A similar comparison has already been performed by the VASP group13 on the basis of the published GAUSSIAN results, and they argued that the most probable explanation for the deviations in Tables 3 and 4 between the PAW and the GTO Table 4. Band Transitions at High-Symmetry Points of the First Brillouin Zone for Silicon, Diamond, and Magnesium Oxidea RIPER PBE

PBE0

Si Γv → Γc Γv → Xc Γv → Lc

2.58 0.78 1.61

6-31G* 3.96 1.98 2.93

C Γv → Γc Γv → Xc Γv → Lc

5.57 4.74 8.47

MgO Γv → Γc Γv → Xc Γv → Lc MgO Γv → Γc Γv → Xc Γv → Lc

VASP HSE06

PBE

PBE0

HSE06

3.33 1.36 2.31

2.57 0.71 1.54

PAW 3.97 1.93 2.88

3.32 1.29 2.24

6-31G* 7.67 6.64 10.78

6.95 5.92 10.04

5.59 4.76 8.46

PAW 7.69 6.66 10.77

6.97 5.91 10.02

5.12 10.35 9.26

Catti 7.68 13.11 11.92

6.93 12.33 11.15

4.75 9.15 7.91

PAW 7.24 11.67 10.38

6.50 10.92 9.64

4.80 9.44 7.94

Catti-mod 7.28 11.93 10.49

6.55 11.18 9.74

4.75 9.15 7.91

PAW 7.24 11.67 10.38

6.50 10.92 9.64

a

We compare our results of PBE, PBE0, and HSE06, obtained with RIPER, to the PAW results of Paier et al.,13,56 determined with VASP. We use the 6-31G* basis sets for Si and C,54 while for MgO we employ the basis set published by Catti et al.,55 referred to as “Catti”, and a corresponding modified one “Catti-mod”, explained in the text. We use the same nonzero-point corrected experimental lattice constants as the VASP group.13 All values in the table are given in eV. J


Article


Figure 6. MgO band structures for primitive cells as calculated with RIPER and QUANTUM ESPRESSO, abbreviated “QE” in the legend. For all calculations MgO is assumed to be in the rocksalt structure at the nonzero-point corrected experimental lattice constant of 4.207 Å.13 Results for RIPER employing (a) the Catti basis set55 and (b) the modified basis set “Catti-mod”, explained in the text. The calculations with QUANTUM ESPRESSO use the PAW method and are identical in both panels.

to speedups of 100−1000 compared to the corresponding Γpoint calculations. The overall performance of calculations involving Fock exchange is dominated by the number of calculated four-center ERIs, which is in turn determined by the size of the BvK supercell as well as the four-center integral neglect threshold. Table 2 indicates the proportionality between CPU time and number of four-center integral batches. This proportionality should not be overemphasized, however, because we count the number of batches, ignoring distance, angular momenta, or contraction of the shells, which all influence the time needed for the calculation of a batch. We have implemented a parallel evaluation of the exchange matrix based on the OpenMP interface58 and thereby perform a self-consistent HF-SCF calculation of the LiH crystal in less than 34 min wall-time, using the large basis set proposed by Paier et al.47 and a single Xeon processor with 20 cores. This timing is obtained by employing a 7 × 7 × 7 k-mesh together with a fourcenter integral neglect threshold of 10−12 a.u. Timings can be further reduced drastically to 15.5 and 5.5 min by using less accurate ERI thresholds of 10−10 a.u. and 10−8 a.u., respectively. In all three calculations, the initial wave function was obtained from a fully converged PBE calculation. Under the same conditions but with the smaller pob-TZVP basis set, the fully converged HF-SCF calculation with a four-center integral neglect threshold of 10−12 a.u. takes less than 2 min wall time.

range up to 0.36% and 8.8% for a0 and B0, they reduce to at most 0.22% and 3.0% for Catti-mod, respectively.

5. PERFORMANCE The computational demand of pure HF calculations or of DFT calculations with a hybrid functional is largely determined by the construction of the Fock exchange matrix K̃ lμν. We therefore present here only timings for the building of K̃ lμν, as obtained on a 2.40 GHz Intel Xeon E5-2680v4 CPU. A performance study of the underlying DFT program RIPER can be found in ref 26. Figure 7 shows CPU timings of the Fock exchange matrix evaluation for bulk LiH using several BvK supercell sizes. In one

Figure 7. Timings for the construction of a single Fock exchange matrix for LiH using the basis set pob-TZVP and lattice constant 4.084 Å (see also Table 1). Black data represent calculations at the Γ-point only, containing 250 to 1458 atoms in the BvK supercell. Red data display simulations, which employ two atoms in the primitive cell and k-meshes with 5 × 5 × 5 up to 27 × 27 × 27 points. This corresponds to BvK supercells consisting of 250 to 39,366 atoms. Solid lines indicate linear (red) and quadratic (black) scaling behavior.

6. CONCLUSIONS We have presented a robust implementation of the periodic Fock exchange using GTOs. We have discussed the divergence of exchange matrix elements K̃ lμν, arising from the artificial periodicity of the density matrix Dlμν = Dl+L μν , which is always introduced in practical calculations by the discretization of the employed set of wavevectors. This discrete set of wavevectors determines in turn the size of the BvK supercell. We have compared two schemes, MIC and TCI, that remedy the divergence for discrete k-meshes and have identified the necessary size of supercells in order to achieve stable SCF calculations and convergence of the total HF energy for selected crystals. While we observe instabilities for small supercells, we recover the same HF energies with MIC and TCI methods for large supercells. The “necessary” size of the supercell or k-mesh for a reliable HF energy is determined by the locality of the

set of calculations, we treat all atoms in the supercells explicitly, working at the Γ-point only. In the other set, we use a primitive cell combined with k-point meshes, which range from 5 × 5 × 5 up to 27 × 27 × 27 k-points. For the Γ-point calculations, a line of slope 2 in the log−log plot indicates the approximate O(N2) scaling of the algorithm. In the k-point-assisted calculation, the scaling is slightly sublinear. Aside from the different scaling behavior, we observe that the k-point simulations show a much lower prefactor. Employing k-point symmetry in our tests leads K



■

ACKNOWLEDGMENTS A.I. and F.P. thank J. C. Klöckner and M. Sierka for stimulating discussions. Both of them also gratefully acknowledge funding from the Carl Zeiss Foundation and from the Collaborative Research Center (SFB) 767 of the German Research Foundation (DFG). A.M.B. thanks the Fonds der Chemischen Industrie for financial support. An important part of the numerical work was carried out on the computational resources of the bwHPC program, namely, the bwUniCluster and the JUSTUS HPC facility, which are supported by the state of Baden-Württemberg and the DFG through Grant INST 40/ 467-1 FUGG.

density matrix and hence depends both on the electronic structure of the studied material and the chosen basis set. For selected insulators and semiconductors, we have demonstrated numerically that the HF total energy converges exponentially with the number of k-points. Our work demonstrates that the MIC approach may be used to calculate the accurate value of the HF energy of the bulk lithium hydride benchmark system, employing a large and flexible basis set. Our results raise doubts that the instability of total energies for increasing basis set quality, reported in ref 19 for two water molecules in a periodic box, is merely based on the MIC scheme. In addition, we have achieved excellent agreement for lattice constants, bulk moduli, and band transitions of Si, diamond, and MgO between values calculated within our MIC scheme and data published in GTO and PAW studies. Although MIC and TCI lead to identical energies for large BvK supercells, stability for medium-size supercells is critical to establish efficient HF and hybrid functional methods in computational practice. We observe that our TCI implementation requires larger supercells compared to our MIC scheme in order to achieve SCF stability and converged HF energies. From our analysis of the density and exchange matrices, we attribute this behavior of TCI to the boundary artifacts in GTO-based four-center integrals near the cutoff radius of TCI. The CRYSTAL code successfully uses overlap and pseudooverlap criteria to truncate the exchange series.59,60 While the pseudo-overlap criterion can be seen as “empirical” and as being of “partially arbitrary nature”, quoting the programmers,59,60 our studies with TCI and MIC show that for large enough BvK supercells the precise truncation scheme becomes irrelevant. We prefer to use MIC and TCI since they basically avoid by our construction the unphysical interaction of periodic images of the crystal orbitals already for intermediate-size BvK supercells. The artifacts in the density matrix elements, seen around the cutoff radius of the Coulomb interaction in TCI-based HF calculations, imply a distortion of the corresponding HF orbitals. The subsequent use in correlated wave function methods, such as Møller−Plesset perturbation theory, needs careful assessment of the impact on accuracy.61 Improvements of the TCI scheme may be achieved by smearing out the cutoff. Progress along these lines was presented in ref 62 for MP2 calculations utilizing a truncated Coulomb kernel, for which the spatial cutoff was varied from smooth to sharp. The efficient evaluation of GTO-based four-center integrals in real space allows the use of large basis sets in combination with fine k-meshes. From our point of view, the next milestone for the widespread use of GTO-based codes in the description of solidstate systems will be the systematic construction of high-quality basis sets. We have shown exemplarily for MgO that a refined GTO basis is able to reproduce both occupied and unoccupied electronic bands with errors that are in the range of present PAW methods.

■

Article

■

REFERENCES

(1) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge University Press: Cambridge, 2004. (2) Martin, R. M.; Reining, L.; Ceperley, D. M. Interacting Electrons: Theory and Computational Approaches; Cambridge University Press: Cambridge, 2016. (3) Becke, A. D. A new mixing of Hartree-Fock and local densityfunctional theories. J. Chem. Phys. 1993, 98, 1372−1377. (4) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys. 1996, 105, 9982−9985. (5) Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A. Density functional theory is straying from the path toward the exact functional. Science 2017, 355, 49−52. (6) Pisani, C.; Maschio, L.; Casassa, S.; Halo, M.; Schütz, M.; Usvyat, D. Periodic local MP2 method for the study of electronic correlation in crystals: Theory and preliminary applications. J. Comput. Chem. 2008, 29, 2113−2124. (7) Marsman, M.; Grüneis, A.; Paier, J.; Kresse, G. Second-order Møller-Plesset perturbation theory applied to extended systems. I. Within the projector-augmented-wave formalism using a plane wave basis set. J. Chem. Phys. 2009, 130, 184103. (8) Grüneis, A. A coupled cluster and Møller-Plesset perturbation theory study of the pressure induced phase transition in the LiH crystal. J. Chem. Phys. 2015, 143, 102817. (9) Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 054115. (10) Grundei, M. M. J.; Burow, A. M. Random Phase Approximation for Periodic Systems Employing Direct Coulomb Lattice Summation. J. Chem. Theory Comput. 2017, 13, 1159−1175. (11) Klimeš, J.; Kaltak, M.; Kresse, G. Predictive GW calculations using plane waves and pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 075125. (12) Pisani, C.; Dovesi, R. Exact-exchange Hartree-Fock calculations for periodic systems. I. Illustration of the method. Int. J. Quantum Chem. 1980, 17, 501−516. (13) Paier, J.; Marsman, M.; Hummer, K.; Kresse, G.; Gerber, I. C.; Á ngyán, J. G. Screened hybrid density functionals applied to solids. J. Chem. Phys. 2006, 124, 154709. (14) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (15) Betzinger, M.; Friedrich, C.; Blügel, S. Hybrid functionals within the all-electron FLAPW method: Implementation and applications of PBE0. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 195117.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Asbjörn M. Burow: 0000-0003-3651-6650 Fabian Pauly: 0000-0001-8017-2379 Notes

The authors declare no competing financial interest. L


Article

Journal of Chemical Theory and Computation (16) Schlipf, M.; Betzinger, M.; Friedrich, C.; Ležaić, M.; Blügel, S. HSE hybrid functional within the FLAPW method and its application to GdN. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 125142. (17) Levchenko, S. V.; Ren, X.; Wieferink, J.; Johanni, R.; Rinke, P.; Blum, V.; Scheffler, M. Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework. Comput. Phys. Commun. 2015, 192, 60−69. (18) Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: Assessment of the Heyd−Scuseria−Ernzerhof screened Coulomb hybrid functional. J. Chem. Phys. 2004, 121, 1187− 1192. (19) Guidon, M.; Hutter, J.; VandeVondele, J. Robust Periodic Hartree-Fock Exchange for Large-Scale Simulations Using Gaussian Basis Sets. J. Chem. Theory Comput. 2009, 5, 3010−3021. (20) Tymczak, C. J.; Weber, V. T.; Schwegler, E.; Challacombe, M. Linear scaling computation of the Fock matrix. VIII. Periodic boundaries for exact exchange at the Γ point. J. Chem. Phys. 2005, 122, 124105. (21) Gygi, F.; Baldereschi, A. Self-consistent Hartree-Fock and screened-exchange calculations in solids: Application to silicon. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 4405−4408. (22) Spencer, J.; Alavi, A. Efficient calculation of the exact exchange energy in periodic systems using a truncated Coulomb potential. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 193110. (23) Sundararaman, R.; Arias, T. A. Regularization of the Coulomb singularity in exact exchange by Wigner-Seitz truncated interactions: Towards chemical accuracy in nontrivial systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 165122. (24) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207− 8215. (25) Izmaylov, A. F.; Scuseria, G. E.; Frisch, M. J. Efficient evaluation of short-range Hartree-Fock exchange in large molecules and periodic systems. J. Chem. Phys. 2006, 125, 104103. (26) Łazarski, R.; Burow, A. M.; Sierka, M. Density Functional Theory for Molecular and Periodic Systems Using Density Fitting and Continuous Fast Multipole Methods. J. Chem. Theory Comput. 2015, 11, 3029−3041. (27) Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; Weigend, F. Turbomole. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 91−100. (28) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; Dover Publications: Mineola, 1996. (29) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133−A1138. (30) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188−5192. (31) Ashcroft, N. W.; Mermin, N. D. Solid State Physics; Harcourt: Orlando, 1976. (32) Causa, M.; Dovesi, R.; Orlando, R.; Pisani, C.; Saunders, V. R. Treatment of the exchange interactions in Hartree-Fock LCAO calculations of periodic systems. J. Phys. Chem. 1988, 92, 909−913. (33) Burow, A. M.; Sierka, M. Linear Scaling Hierarchical Integration Scheme for the Exchange-Correlation Term in Molecular and Periodic Systems. J. Chem. Theory Comput. 2011, 7, 3097−3104. (34) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (35) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 2006, 125, 224106. (36) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)]. J. Chem. Phys. 2006, 124, 219906. (37) Pisani, C.; Aprà, E.; Causà, M. Density matrix of crystalline systems. I. Long-range behavior and related computational problems. Int. J. Quantum Chem. 1990, 38, 395−417. (38) Carrier, P.; Rohra, S.; Görling, A. General treatment of the singularities in Hartree-Fock and exact-exchange Kohn-Sham methods for solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 205126.

(39) Stolarczyk, L. Z.; Piela, L. Direct Calculation of Lattice Sums. A Method to Account for the Crystal Field Effects. Int. J. Quantum Chem. 1982, 22, 911−927. (40) Kohn, W. Analytic Properties of Bloch Waves and Wannier Functions. Phys. Rev. 1959, 115, 809−821. (41) Kohn, W. Density functional theory for systems of very many atoms. Int. J. Quantum Chem. 1995, 56, 229−232. (42) Goedecker, S. Decay properties of the finite-temperature density matrix in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3501−3502. (43) Ismail-Beigi, S.; Arias, T. A. Locality of the Density Matrix in Metals, Semiconductors, and Insulators. Phys. Rev. Lett. 1999, 82, 2127−2130. (44) Ahlrichs, R. A simple algebraic derivation of the Obara−Saika scheme for general two-electron interaction potentials. Phys. Chem. Chem. Phys. 2006, 8, 3072−3077. (45) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular ElectronicStructure Theory; John Wiley & Sons, Inc.: Chichester, 2013. (46) Peintinger, M. F.; Oliveira, D. V.; Bredow, T. Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations. J. Comput. Chem. 2013, 34, 451−459. (47) Paier, J.; Diaconu, C. V.; Scuseria, G. E.; Guidon, M.; VandeVondele, J.; Hutter, J. Accurate Hartree-Fock energy of extended systems using large Gaussian basis sets. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 174114. (48) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (49) Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chem. Phys. Lett. 1980, 73, 393−398. (50) Häser, M.; Ahlrichs, R. Improvements on the direct SCF method. J. Comput. Chem. 1989, 10, 104−111. (51) Civalleri, B.; Orlando, R.; Zicovich-Wilson, C. M.; Roetti, C.; Saunders, V. R.; Pisani, C.; Dovesi, R. Comment on “Accurate HartreeFock energy of extended systems using large Gaussian basis sets”. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 106101. (52) Murnaghan, F. D. The compressibility of media under extreme pressure. Proc. Natl. Acad. Sci. U. S. A. 1944, 30, 244−247. (53) Grüneis, A.; Marsman, M.; Kresse, G. Second-order Møller− Plesset perturbation theory applied to extended systems. II. Structural and energetic properties. J. Chem. Phys. 2010, 133, 074107. (54) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. J. Chem. Phys. 1982, 77, 3654−3665. (55) Catti, M.; Valerio, G.; Dovesi, R.; Causà, M. Quantummechanical calculation of the solid-state equilibrium MgO+αAl2O3⇄MgAl2O4 (spinel) versus pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14179−14187. (56) Paier, J.; Marsman, M.; Hummer, K.; Kresse, G.; Gerber, I. C.; Á ngyán, J. G. Erratum: “Screened hybrid density functionals applied to solids” [J. Chem. Phys. 124, 154709 (2006)]. J. Chem. Phys. 2006, 125, 249901. (57) We used the projectors Mg.pbe-nsp-bpaw.UPF and O.pbekjpaw.UPF from http://www.quantum-espresso.org (accessed August 2018). (58) OpenMP Architecture Review Board, OpenMP Application Program Interface, Version 4.0, July 2013. http://www.openmp.org (accessed August 2018). (59) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock Ab Initio Treatment of Crystalline Systems; Lecture Notes in Chemistry; Springer: Berlin, 1988; Vol. 48. (60) Roetti, C. In Quantum-Mechanical Ab-initio Calculation of the Properties of Crystalline Materials; Lecture Notes in Chemistry; Pisani, C., Ed.; Springer: Berlin, 1996; Vol. 67. (61) Del Ben, M.; Hutter, J.; VandeVondele, J. Electron Correlation in the Condensed Phase from a Resolution of Identity Approach Based on the Gaussian and Plane Waves Scheme. J. Chem. Theory Comput. 2013, 9, 2654−2671. M


Article

Journal of Chemical Theory and Computation (62) Dutoi, A. D.; Head-Gordon, M. A Study of the Effect of Attenuation Curvature on Molecular Correlation Energies by Introducing an Explicit Cutoff Radius into Two-Electron Integrals. J. Phys. Chem. A 2008, 112, 2110−2119.

N


Robust Periodic Fock Exchange with Atom-Centered Gaussian Basis

Recommend Documents