Analytical Approach to Phonon Calculations in the SCC-DFTB Framework
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug Analytical Approach to Phonon Calculations in the SCC-DFTBFramework
Vladimir Baˇci´c, Thomas Heine , , and Agnieszka Kuc Department of Physics and Earth Sciences, Jacobs University Bremen,Campus Ring 1, 28759 Bremen, Germany Helmholtz-Zentrum Dresden-Rossendorf, Abteilung Ressourcen¨okologie, Forschungsstelle Leipzig,Permoserstr. 15, 04318 Leipzig, Germany Theoretical Chemistry, TU Dresden, Mommsenstr. 13, 01062 Dresden, Germany Department of Chemistry, Yonsei University,50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Korea
Abstract
Detailed derivation of the analytical, reciprocal-space approach of Hessian calculation within the self-consistent-charge density functional based tight-binding framework (SCC-DFTB) is presented. This ap-proach provides an accurate and efficient way for obtaining the SCC-DFTB Hessian of periodic systems. Itssuperiority with respect to the traditional numerical force differentiation method is demonstrated for dopedgraphene, graphene nanoribbons, boron-nitride nanotubes, bulk zinc-oxide and other systems.
In the past two decades, density-functional based tight-binding method (DFTB) has become a relativelypopular tool for quantum mechanical simulations of large systems, otherwise computationally too demandingfor the standard density-functional theory (DFT) or other ab-initio methods. Introduced in the mid-90s asan approximation to DFT, DFTB has been subjected to ongoing extensions and improvements, the so-calledself-consistent charge (SCC) DFTB being perhaps the most important one. In spite of its shortcomings, DFTB(along with its extensions) has shown to perform reasonably well for a variety of systems, sometimes even with anaccuracy comparable to that of its first-principles counterpart, with only a fraction of the computational cost. Having this in mind, and motivated by the need of having an efficient and reasonably accurate way for studyingthe vibrational properties of large systems, Witek et al., and subsequently Nishimoto and Irle, developed ananalytical method for obtaining the Hessian (i.e., matrix of geometrical second derivatives, needed in calculationof vibrational frequencies), within the SCC-DFTB framework. Although considerably more efficient than thetraditional numerical force differentiation method of Hessian calculation, the application of the approach fromrefs. and to periodic systems requires using supercells, which results in lower overall accuracy and highercomputational costs.In this paper, we present an analytical and supercell-free method for calculating the SCC-DFTB Hessian ofperiodic systems. The underlying approach is based on the direct evaluation of the discrete Fourier transform ofthe Hessian (rather than the Hessian itself), which can be achieved by taking the second derivative of the totalenergy with respect to collective and phase-modulated atomic displacements. This reciprocal-space method ofphonon calculations is not new, in fact, it has been implemented in the plane-wave and muffin-tin based DFTcodes long time ago, and is also known as the density-functional perturbation theory (DFPT) or linear-response(LR) theory. The objective of this work is to develop the SCC-DFTB analogue of such an approach to phononcalculations.This paper is organized as follows: Section 2.1 provides a short review of the standard DFTB equations,along with the notation used throughout the rest of the paper. In Section 2.2, we discuss the main ideabehind the reciprocal-space approach to Hessian calculation. In Section 2.3, we present a derivation of theexpression for evaluating the Fourier-transformed Hessian within SCC-DFTB. For the sake of clarity, only themost important steps are shown here, while complete and extensive mathematical details can be found in theSupplementary Material. Comparison of the reciprocal-space and the numerical force differentiation method ofphonon calculation, along with the underlying discussion, is given in Section 3, while some concluding remarksare given in Section 4.The theoretical formalism developed in this work has been implemented in a locally modified version of DFTBprogram in Amsterdam Modeling Suite, version 2018.1 Theory
We begin by giving a brief overview of the SCC-DFTB framework, without going too deep into details. Amore thorough discussion on this topic, as well as the derivations of the equations presented in this Section,can be found in the literature.
In the SCC-DFTB approximation, the Kohn-Sham (KS) energy functionalfor periodic systems can be written as: E (cid:2) { c k ,n } (cid:3) = X k ,n f k ,n c † k ,n H k c k ,n + 12 X I,J X R γ IJ ( R ) ∆z I ∆z J + 12 X I,J X R V rplIJ ( R ) (1)In the first term, k is a vector in the Brillouin zone, n is the band index, c k ,n is a column-vector of the orbitalcoefficients, H k is the parametrized Hamiltonian matrix (in the Bloch basis) and f k ,n is the electron occupationfunction. In the last two terms, the sum runs over all atom pairs ( I, J ) and all lattice vectors R . In the secondterm, the function γ IJ describes the Coulomb interaction of atomic charge fluctuations ∆z . Finally, V rpl is theso-called repulsion potential, parametrized as a short-ranged isotropic force-field. The three terms in eq. (1)are called band-structure ( E BS ), charge-fluctuation ( E CF ) and repulsion energy ( E rpl ), respectively. We notethat H k , γ IJ and V rpl depend only on the pre-calculated parameter set and the geometry of the system, butnot on the orbital coefficients.The charge fluctuations in (1) are most commonly calculated using the Mulliken population analysis. According to it, the charge fluctuation on atom I is given by: ∆z I = 12 X k ,n f k ,n X a ∈ I X b (cid:0) c a ∗ k ,n c b k ,n S ab k + c b ∗ k ,n c a k ,n S ba k (cid:1) − z I (2)where z I is the valence charge of the corresponding neutral atom, S ab k is the overlap matrix element betweenatomic basis functions a and b (also in the Bloch basis). To avoid explicitly writing the double sum over thebasis functions, it is convenient to introduce the atom projection matrix P I and the charge projection matrix Z I, k , defined by: P abI ≡ ( δ a,b , if basis functions a and b belong to atom I Z I, k ≡
12 ( P I S k + S k P I ) (4)Using (3) and (4), the expression for the Mulliken charge fluctuations (2) can be compactly written as: ∆z I = 12 X k ,n f k ,n c † k ,n ( P I S k + P I S k ) c k ,n − z I = X k ,n f k ,n c † k ,n Z I, k c k ,n − z I (5)According to the variational principle, the SCC-DFTB ground state energy is obtained by minimizing (1)with respect to orbital coefficients, under the orthonormalization constraints: c † k ,m S k c k ,n = δ m,n . This leadsto a system of generalized eigenvalue equations for c k ,n : h H k + X I V I Z I, k i c k ,n = ε k ,n S k c k ,n (6)where the eigenvalues ε k ,n are the single-particle band-structure energies, and V I ≡ X J X R γ IJ ( R ) ∆z J (7)is the electrostatic potential on atom I due to charge fluctuations on all atoms. The term in the square bracketsof (6) can be regarded as the total Hamiltonian matrix H k . As H k includes the term with the charge fluctu-ations, which depend on the orbital coefficients, equations (5) and (6) have to be solved self-consistently, justlike the standard DFT KS equations. The SCC-DFTB ground-state energy is then given by evaluating theexpression (1) with the self-consistent orbital coefficients and charge fluctuations.2 .2 Energy derivatives and vibrational properties In the Born-Oppenheimer approximation, the frequencies and modes of phonons with wavevector q aregiven as the eigenvalues and eigenvectors of the so-called dynamical matrix D q , defined as:D ( B,β ; A,α ) q ≡ √ M B M A X R e i qR Φ ( B,β ; A,α ) R (8)where capital and Greek indexes denote atoms and Cartesian displacements, respectively, M is the atomic massand: Φ ( B,β ; A,α ) R ≡ ∂ E∂u βB, R ∂u αA, (9)where u µX, R denotes the µ -th Cartesian component of the position of atom X belonging to lattice point R , isknown as the interatomic force constant matrix, Hessian matrix or simply Hessian. Although Φ R is formallydefined for all points of the Bravais lattice { R } of the system, in practice, it has non-negligible values only onsome finite subset of { R } , i.e., on a supercell of the underlying system. Since the dynamical matrix is justthe discrete Fourier transform of the Hessian, weighted by the inverse square root of atomic mass products,calculating the Hessian poses the main challenge in the study of vibrational properties and related phenomena.The conceptually easiest approach to this problem consists of numerical evaluation of the first-order forcederivatives with respect to atomic displacements. Although simple, this technique of Hessian calculation canbe quite slow and inefficient, since it requires doing a number of force calculations on a supercell on which theHessian is non-negligible, thereby typically resulting in much higher computational costs compared to calcula-tions on the corresponding primitive unit cell. Numerical instabilities associated with numerical evaluation ofderivatives can also present a more severe issue in this case. Even if an analytical expression for evaluating theforce derivatives is available, the problem of using supercells is still present.An alternative approach to phonon calculation is based on evaluating the discrete Fourier transform of theHessian directly, namely by using the following identity: e Φ ( B,β ; A,α ) q = ˜ ∂ B,β - q ˜ ∂ A,α q E (10)where ˜ ∂ X,µ q ≡ X R e i qR ∂∂u µX, R (11)is the discrete Fourier-transform of the position derivative operator. Since (10) is valid for any q -point ofthe reciprocal space, the entire phonon spectrum can be obtained without calculating Φ R at all. In practice,however, it is generally much more efficient to evaluate e Φ q on a regular grid of the Brillouin zone (oftenreferred to as the q-grid), apply the inverse Fourier transformation to get Φ R , and then use (8) to compute thedynamical matrix at arbitrary q -point. From the properties of the discrete Fourier transform, it follows thatthe density of the q-grid determines the size of the supercell on which Φ R is defined. Whether there is anactual advantage to this (i.e., reciprocal-space) approach to Hessian calculation, depends on the way the totalenergy is calculated. In other words, for this approach to be useful for a given energy calculation method, onemust be able to efficiently evaluate the right-hand side (RHS) of (10) within that particular method. In thefollowing Section, we shall see that within SCC-DFTB, this can in fact be done analytically and without usingsupercells , just like in plane-wave and mixed-basis DFT formalisms.Before moving on, we briefly focus on the Fourier-transformed atomic position derivative operator (11). If F R is any atomic-position-dependent quantity defined on all lattice points of a periodic system (e.g., orbitalcoefficients, charge fluctuations etc.), then acting with ˜ ∂ q on it results in a phase-modulated quantity:˜ ∂ A,µ q F R = e i qR X R ′ e i qR ′ ∂F R ∂u µA, R + R ′ ≡ e i qR F ( A, q ) R (12)But no such modulation is present when acting on F R with two Fourier-transformed derivative operators withopposite wave-vectors: ˜ ∂ B,ν - q ˜ ∂ A,µ q F R = F ( B, - q ; A, q ) R (13)On the far RHS of (12) and (13), the indexes of Cartesian displacements in the superscript have been droppedfor clarity. Unless specified otherwise, such notation shall be used from now on, i.e., any quantity with asuperscript containing the atom index and wave-vector pair(s) shall be assumed to be the phase-independentpart of the quantity obtained by applying the Fourier-transformed derivative operator(s).3 .3 Second derivatives of SCC-DFTB energy In this Section, we derive the expression for evaluating the RHS of (10), for the case where E is the SCC-DFTB total energy. Following refs. and, we do this by applying the variational principle to the secondderivative of the DFTB KS energy functional (1):˜ ∂ B,β - q ˜ ∂ A,α q E = min h ˜ ∂ B,β - q ˜ ∂ A,α q E i (14)This way, the expressions for the second derivatives of each component of E (band-structure, charge-fluctuationand repulsion) can be derived separately, so this is how we proceed. To keep the discussion as clear and assimple as possible, only systems with a finite band-gap and integer electron occupations (i.e., insulators andsemiconductors at zero electron temperature) shall be considered here, while the corresponding equations forthe general case are given in the Supplementary Material. A Second derivatives of E BS The second derivative of the band-structure part of E reads:˜ ∂ B - q ˜ ∂ A q E BS = X k ,n f k ,n (cid:20) c ( B, - q ; A, q ) † k ,n H k c k ,n + c † k ,n H B, - q ; A, q ) k c k ,n + c † k ,n H k c ( B, - q ; A, q ) k ,n + c ( B, q ) † k ,n H A, q ) k c k ,n + c ( B, q ) † k ,n H k + q c ( A, q ) k ,n + c † k ,n H B, q ) † k c ( A, q ) k ,n + c ( A, - q ) † k ,n H B, - q ) k c k ,n + c ( A, - q ) † k ,n H k − q c ( B, - q ) k ,n + c † k ,n H A, - q ) † k c ( B, - q ) k ,n (cid:21) (15)The matrix elements of H A, q ) k and H B, - q ; A, q ) k are given by: h H A, q ) k i ab = X R ∇ H ab ( R + u ab ) e i kR (cid:16) − δ Aa + e i qR δ Ab (cid:17) (16a) h H B, - q ; A, q ) k i ab = X R ∇ H ab ( R + u ab ) e i kR (cid:16) δ A,B (cid:0) δ Aa + δ Ab (cid:1) − e − i qR δ Aa δ Bb − e i qR δ Ab δ Ba (cid:17) (16b)where u ab is the bond vector between atoms to which basis functions a and b belong, and δ Xx =1 if basis function x belongs to atom X , zero otherwise. A completely equivalent expression holds for the matrix elements of S ( A, q ) k and S ( B, - q ; A, q ) k as well. It is also worth mentioning that the RHS in (16) can in fact be evaluated analytically(see appendix A), which contributes to the overall efficiency and accuracy of the reciprocal-space approach. B Second derivatives of E CF The second derivative of the charge-fluctuation part of E reads:˜ ∂ B - q ˜ ∂ A q E CF = X I,J h ∆z ( A, q ) I γ ( B, - q ) IJ ∆z J + ∆z ( B, - q ) I γ ( A, q ) IJ ∆z J + ∆z ( B, - q ) I e γ IJ ( q ) ∆z ( A, q ) J + ∆z I e γ IJ ( ) ∆z ( B, - q ; A, q ) J i + ˜ ∂ B - q ˜ ∂ A q E CF [ { ∆z } ] (17)Here, e γ IJ ( q ) is the phase-modulated lattice sum of γ IJ ( R ): e γ IJ ( q ) ≡ X R γ IJ ( R ) e i qR (18) γ ( A, q ) IJ is defined as: γ ( A, q ) IJ ≡ X R (cid:16) − δ J,A ∇ γ IA ( R ) e i qR + δ I,A ∇ γ AJ ( R ) (cid:17) (19)and: ˜ ∂ B - q ˜ ∂ A q E CF [ { ∆z } ] = ∆z A X R (cid:16) − ∆z B ∇ γ BA ( R ) e i qR + δ A,B X I ∆z I ∇ γ AI ( R ) (cid:17) (20)is the second derivative of the charge fluctuation interaction with respect to atom positions only (i.e., keepingthe charge fluctuations constant). The problem in (18)-(20) and (7) is that γ is a long-ranged function (itdecays as slowly as the Coulomb potential ), which makes the underlying lattice summations only conditionallyconvergent. To overcome this issue, we use the well-known Ewald summation technique, more details on thiscan be found in Appendix B. 4hase-independent part of the charge-fluctuation first derivative can be written in terms of orbital coefficientsand their first derivatives as: ∆z ( A, q ) I = 12 X k ,n f k ,n h c † k ,n (cid:0) P I S k + q + S k P I (cid:1) c ( A, q ) k ,n + c ( A, - q ) † k ,n (cid:0) P I S k + S k − q P I (cid:1) c k ,n i + ∆z ( A, q ) I (21)with the last term defined as: ∆z ( A, q ) I ≡ X k ,n f k ,n c † k ,n (cid:16) P I S ( A, q ) k + S ( A, - q ) † k P I (cid:17) c k ,n (22)Finally, the second derivative of the charge-fluctuations is given by: ∆z ( B, - q ; A, q ) I = 12 X k ,n f k ,n h c ( B, - q ; A, q ) † k ,n Z I, k c k ,n + c † k ,n Z I, k c ( B, - q ; A, q ) k ,n + c † k ,n Z ( B, - q ; A, q ) I, k c k ,n + c ( B, q ) † k ,n Z ( A, q ) I, k c k ,n + c ( B, q ) † k ,n Z I, k + q c ( A, q ) k ,n + c † k ,n Z ( B, q ) † I, k c ( A, q ) k ,n + c ( A, - q ) † k ,n Z ( B, - q ) I, k c k ,n + c ( A, - q ) † k ,n Z I, k − q c ( B, - q ) k ,n + c † k ,n Z ( A, - q ) † I, k c ( B, - q ) k ,n i (23)where Z ( A, q ) I, k and Z ( B, - q ; A, q ) I, k are given exactly as in (4), with S k replaced by S ( A, q ) k and S ( B, - q ; A, q ) k , respectively. C Second derivatives of E rpl The second derivative of the repulsion energy is given by:˜ ∂ B - q ˜ ∂ A q E rpl = X R (cid:16) − ∇ V rplBA ( R ) e i qR + δ A,B X I ∇ V rplAI ( R ) (cid:17) (24)In contrast to γ , the DFTB repulsion potential is generally a short-ranged function which decays rapidly withdistance, so evaluating the lattice on the RHS of (24) is straightforward. Just like the repulsion energy, therepulsion contribution to the Hessian can be obtained independently from the band-structure and charge-fluctuation contributions. Final expressions
So far, we have derived expressions for the second derivatives of all components of E . To proceed further,the second derivatives of the orthonormalization constraints must be taken into account as well: c ( B, - q ; A, q ) † k ,n S k c k ,n + c † k ,n S ( B, - q ; A, q ) k c k ,n + c † k ,n S k c ( B, - q ; A, q ) k ,n + (cid:16) c ( B, q ) † k ,n S ( A, q ) k c k ,n + c ( B, q ) † k ,n S k + q c ( A, q ) k ,n + c † k ,n S ( B, q ) † k c ( A, q ) k ,n (cid:17) + (cid:16) A ↔ B q ↔ − q (cid:17) = 0 (25)The left-hand side (LHS) of this equation is fully analogous to the term in the square brackets of (15), withthe Hamiltonian replaced by the overlap matrix.Now, adding (15), (17) and (24), while making use of (18)-(23) and (25), the second derivative of the SCC-DFTBenergy functional becomes:˜ ∂ B - q ˜ ∂ A q E = X k ,n f k ,n (cid:20) c ( B, q ) † k ,n (cid:16) H A, q ) k + X I V I Z ( A, q ) I, k − ε k ,n S ( A, q ) k (cid:17) c k ,n + c † k ,n (cid:16) H A, - q ) k + X I V I Z ( A, - q ) I, k − ε k ,n S ( A, - q ) k (cid:17) c ( B, - q ) k ,n + c ( B, q ) † k ,n (cid:16) H k + q − ε k ,n S k + q (cid:17) c ( A, q ) k ,n (cid:21) + (cid:20) A ↔ B q ↔ − q (cid:21) + X I,J (cid:18) γ ( B, - q ) IJ ∆z ( A, q ) I ∆z J + γ ( A, q ) IJ ∆z ( B, - q ) I ∆z J + e γ IJ ( q ) ∆z ( B, - q ) I ∆z ( A, q ) J (cid:19) + ˜ ∂ B - q ˜ ∂ A q E BS [ { c k ,n } ] + ˜ ∂ B - q ˜ ∂ A q E CF [ { ∆z } ] + ˜ ∂ B - q ˜ ∂ A q E rpl (26)where: ˜ ∂ B - q ˜ ∂ A q E BS [ { c k ,n } ] ≡ X k ,n f k ,n c † k ,n (cid:16) H B, - q ; A, q ) k + X I V I Z ( B, - q ; A, q ) I, k − ε k ,n S ( B, - q ; A, q ) k (cid:17) c k ,n (27)5ust like ˜ ∂ B - q ˜ ∂ A q E CF [ { ∆z } ], depends only on the orbital coefficients (but not on their derivatives!) and thegeometry of the system.We see that the expression for ˜ ∂ B - q ˜ ∂ A q E , as given by (26), contains no second derivatives of either the orbitalcoefficients or the charge fluctuations. Furthermore, it is variational with respect to the first derivatives of theorbital coefficients and, provided that the orbital coefficients on its RHS minimize E , its minimum correspondsto the true value of the second derivative of the SCC-DFTB energy. Varying (26) results in the following equation for c ( A, ± q ) k ,n : − (cid:16) H k ± q − ε k ,n S k ± q (cid:17) c ( A, ± q ) k ,n = (cid:16) H ( A, ± q ) k − ε k ,n S ( A, ± q ) k (cid:17) c k ,n (28)and an equivalent one for c ( B, ± q ) k ,n . H ( A, ± q ) k is the matrix of the Hamiltonian total derivative, given by: H ( A, ± q ) k ≡ H A, ± q ) k + X I h V I Z ( A, ± q ) I, k + V ( A, ± q ) I (cid:0) P I S k + S k ± q P I (cid:1)i (29)with: V ( A, ± q ) I ≡ X J (cid:16) γ ( A, ± q ) IJ ∆z J + e γ IJ ( ± q ) ∆z ( A, ± q ) J (cid:17) (30)being the total derivative of the electrostatic potential. In the literature, (28) is also known as the Sternheimerequation. So the problem of calculating ˜ ∂ B - q ˜ ∂ A q E effectively reduces to the problem of determining the orbital coefficientderivatives. These can be expressed as linear combinations of the orbital coefficients: c ( X, ± q ) k ,n = X m U ( X, ± q ) k ( n,m ) c k ± q ,m (31)where U ( X, ± q ) k is a square matrix (with the number of rows and columns equal to the number of atomic basisfunctions) to be determined. From (28), it immediately follows that the entries of U ( X, ± q ) k matrix, which referto the non-degenerate pairs of states at k and k ± q points, are given by: U ( X, ± q ) k ( n,m ) = c † k + q ,m (cid:16) H ( A, q ) k − ε k ,n S ( A, q ) k (cid:17) c k ,n ε k ,n − ε k ± q ,m (32)while for determining all other entries of U ( X, ± q ) k , the following relation can be used (see S1): U ( X, ∓ q ) ∗ k ± q ( m,n ) + U ( X, ± q ) k ( n,m ) + c † k + q ,m S ( A, q ) k c k ,n = 0 (33)Inserting (31) to (21) and making use of (32) and (33), the following expression for charge-fluctuation derivativesis obtained: ∆z ( A, q ) I = X k X n ∈V X m ∈C f k ,n M ( A, q ) k ( m,n ) ε k ,n − ε k + q ,m − X m ∈V f k + q ,m O ( A, q ) k ( m,n ) ! c † k ,n (cid:0) P I S k + q + S k P I (cid:1) c k + q ,m + ∆z ( A, q ) I (34)Here, V and C refer to the sets of valence and conduction (i.e., occupied and empty) states, respectively,while M ( A, q ) k ( m,n ) is the (generalized) electron-phonon matrix element and O ( A, q ) k ( m,n ) is the overlap derivative matrixelement: M ( A, q ) k ( m,n ) ≡ c † k + q ,m (cid:16) H ( A, q ) k − ε k ,n S ( A, q ) k (cid:17) c k ,n O ( A, q ) k ( m,n ) ≡ c † k + q ,m S ( A, q ) k c k ,n (35)It is easy to see that electron-phonon matrix elements and charge-fluctuation derivatives depend on each other,so they must be calculated self-consistently, much like the charge-fluctuations and the orbital coefficients. Atlast, combining (26) with (31)-(33), we arrive at the expression for the second derivatives of the SCC-DFTBenergy:˜ ∂ B - q ˜ ∂ A q E = X k X n ∈V (cid:20) X m ∈C f k ,n M ( A, q ) k ( m,n ) M ( B, q ) ∗ k ( m,n ) ε k ,n − ε k + q ,m − X m ∈V f k + q ,m (cid:16) M ( A, q ) k ( m,n ) O ( B, q ) ∗ k ( m,n ) + M ( B, q ) ∗ k ( m,n ) O ( A, q ) k ( m,n ) (cid:17) (cid:21) − X I,J e γ IJ ( q ) ∆z ( B, q ) ∗ I ∆z ( A, q ) J + X I h ∆z ( B, q ) ∗ I V ( A, q ) I + ∆z ( A, q ) ∗ I V ( B, q ) I i + ˜ ∂ B - q ˜ ∂ A q E BS [ { c k ,n } ] + ˜ ∂ B - q ˜ ∂ A q E CF [ { ∆z } ] + ˜ ∂ B - q ˜ ∂ A q E rpl (36)6his expression, along with (34), is the main result of this paper. Since all of the derivatives appearing in thisSection can be evaluated analytically, the entire SCC-DFTB reciprocal-space approach to Hessian calculationcan be considered analytical. Although all expressions here are derived for periodic systems, they are alsovalid for non-periodic systems as well. This can be seen by taking the limit of infinitely large unit cells, thusrestricting all real-space summations to a single unit cell and all k and q -points to the Γ-point. In that case,reciprocal-space summations can be omitted, Fourier-transformed derivatives reduce to ordinary derivativesand our entire formulation becomes equivalent to the one developed by Witek et al. In closing of this section,we once again point out that both (36) and (34) are only valid for systems with a finite band gap and forvanishing electron temperature, while the corresponding expressions for the case of arbitrary temperature canbe found in the Supplementary Material ((S2.9) and (S2.11)).7
Test Calculations
To test the performance of the reciprocal-space (analytical) approach to Hessian calculation and compareit to the traditional numerical force differentiation method, we used both approaches to compute the Hessiansfor a variety of systems of all dimensions. In order to check how much the numerical and analytical resultsdiffer, we calculated the so-called root mean squared relative percentage difference of the Hessians resultingfrom the two approaches: ∆ A , N ≡ vuut N at ) N X R X I,J Φ ( I ; J ) A , R Φ ( I ; J ) N , R − ! · Φ A and Φ N are analytically and numerically obtained Hessians, respectively, N at is the number of atomsin the system and N is the number neighbouring unit cells within the supercell on which the Hessians aredefined. In the analytical approach, this supercell is effectively determined by the q-grid used in the underlyingcalculation (as already mentioned in Section 2.2), hence (37) (and comparing Hessians in general) only makessense if the q-grid parameters used in the calculation of Φ A are equal to the supercell parameters in thecalculation of Φ N . To make the comparison of both methods as consistent as possible, the k-grid parametersused in the reciprocal-space integration were made inversely proportional to the (super)cell size. For example,if the Brillouin zone of some system was sampled with a 12 × × × × × × × × × × Results and Discussion
Table 1 shows examples of ∆ A , N for cubic boron-nitride (zinc-blende phase), where the numerical Hessianswere calculated with a different number of steps and step sizes. We see that ∆ A , N always drops when the stepsize is decreased, whereas such a clear trend is not present when increasing the number of displacement steps.The latter behavior can be attributed to the anharmonic effects, which are not captured by the analyticalapproach at this level of theory, but can always appear in the numerical approach for sufficiently large stepsizes. In any case, it is clear that in the limit of small displacement step sizes, the numerically obtained Hessiansconverge to the analytically obtained one. Similar behavior of ∆ A , N is obtained for all other systems consideredhere, which confirms the accuracy of the analytic approach.Table 1: Root mean square relative percentage difference between analytically and numerically obtained Hes-sians (∆ A , N , see eq. (37)) for zinc-blende BN . Number of steps and step sizes refer to the parameters used inthe numerical calculations. total ).When considering the overall efficiency of any quantum-chemical computational method, it is also importantto take the aspects of parallelizability and symmetry into account. Since our implementation of the analyticalmethod is currently serial and cannot make use of symmetry, all calculations (both numerical and analytical)were carried out on a single CPU core, while disregarding all possible symmetries of the investigated systems.But even if the same calculations had been performed using a parallel implementation that does supportsymmetry, we expect that similar timing ratios would have been obtained, as both approaches are in principleparallelizable and can exploit symmetry to the same extent. In the numerical case, for example, each forceconstant can be evaluated independently, while the number of independent force constants is determined bythe symmetry of the system. Likewise, in the analytical case each element of the e Φ q matrix can be evaluatedindependently, while the number of independent q-grid vectors is also determined by the symmetry of the8ystem. Finally, fig. 1 shows examples of phonon spectra and phonon density of states of some systems from table 2,calculated with SCC-DFTB.Table 2: Hessian calculation timings (in seconds) for the analytical ( t A ) and numerical ( t N ) approaches. N at is the number of atoms in the primitive unit cell of the given system and references refer to the parameter setused in the calculations. system Ref. N at supercell(q-grid) SCC-DFTB DFTB0 t N t A t N /t A t N t A t N /t A fullerene
60 N/A 71.25 7.87 9.05 34.85 3.78 9.22graphene nanoribbon
52 2 3129.82 945.69 3.31 623.13 373.32 1.673 6781.11 1697.48 3.99 1438.13 567.75 2.534 13602.81 2268.69 5.99 3118.92 748.01 4.17BN (18,18)-nanotube
72 2 5908.53 917.00 6.44 1774.76 478.60 3.713 13505.67 2070.13 6.52 4358.30 719.79 6.054 23206.95 2353.30 9.86 7383.50 968.84 7.62BeCl trilayer
12 2 × × ×
32 1 × × × × × × × × × × × × × × × × × × × × × M K ΓA L H A ω ( T H z ) SiC (4H)
Γ M K Γ ω ( T H z ) Si-doped grahpene
Γ X ω ( T H z ) BN (18,18)-nanotube
Figure 1: Phonon dispersions and DOS of selected systems
In summary, we have successfully derived the reciprocal-space approach of Hessian calculation within theSCC-DFTB framework. This approach allows for the Hessian of periodic systems to be obtained accurately andwithout doing any calculations on supercells, while also providing some information about electron-phonon in-teractions, which can be of great importance in the study of transport phenomena. The formulation presentedin this paper effectively generalizes previous analytical methods of SCC-DFTB Hessian calculation, which areless suited for periodic systems. Its efficiency has been demonstrated by performing test calculations on varioussystems of all dimensions, where it showed to be significantly faster than the numerical force-differentiationmethod, especially for large systems.We believe that further research on DFTB-based phonon calculation methods can open a gateway to anefficient ab-initio description of phonon-related properties (such as electrical or thermal conductivity, Ramanspectra, band-gap renormalization, superconductivity) of large systems, otherwise computationally too demand-ing for ordinary DFT methods, and perhaps even lead to more extensive DFTB parameters develpoment forsolid-state applications. This paper can be regarded as a first step in that direction. In our future work, weplan various extensions of the theoretical formulation presented here (such as to DFTB3 framework and addingsupport for spin-polarization), as well as its application to problems of interest.
Acknowledgements
We thank Pier Philipsen, Mirko Franchini, Robert R¨uger, Augusto Oliveira, M. S. Ramzan and Stan vanGisbergen for useful discussions and technical support. The financial support by the Deutsche Forschungsge-meinschaft (GRK 2247/1 (QM3)) is greatly appreciated.
Supporting Information
The data that supports the findings of this study are available within the article [and its supplementarymaterial].
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonablerequest.
Appendix A Derivatives of the DFTB matrix elements
According to the Slater-Koster transformation rules, the two-center integrals I ab , between basis functions a and b , located on atoms separated by vector r , can be written as: I ab ( r ) = X τ R τab ( r ) A τl ( ab ) (ˆ r ) (A.1)where τ is the index of the Slater-Koster integral, r ≡ k r k , and ˆ r ≡ r /r . In the DFTB formalism, the radialfunctions R τab ( r ) are generally given on a numerical grid. However, they can always be cast to an analyticalform, for example, by spline interpolation. Unlike R τab ( r ), the functions A τl ( ab ) (ˆ r ) depend only on the angular10omenta of a and b basis functions and not on their radial shape. Since they are given in a purely analyticalform, their derivatives can be obtained easily. For example, if a is an s - and b a p -type function, we have: A σsp i ( ˆr ) = r i r (A.2a) ∂ j A σsp i ( ˆr ) = δ ij r − r i r j r (A.2b) ∂ k ∂ j A σsp i ( ˆr ) = − r ( δ ij r k + δ ik r j + δ jk r i ) + 3 r i r j r k r (A.2c)where r i is the i-th component of r and ∂ i ≡ ∂/∂r i . Similar expressions can be derived for all other combinationsof angular momenta.Finally, the first and second derivative of I ab can be written as: ∂ i I ab ( r ) = X τ (cid:20) ∂ r R τab r i r A τl ( ab ) + R τab ∂ i A τl ( ab ) (cid:21) (A.3a) ∂ j ∂ i I ab ( r ) = X τ (cid:20) (cid:18) ∂ r R τab − ∂ r R τab r (cid:19) r i r j r A τl ( ab ) + R τab ∂ i ∂ j A τl ( ab ) + ∂ r R τab r (cid:16) δ ij A τl ( ab ) + r i ∂ j A τl ( ab ) + r j ∂ i A τl ( ab ) (cid:17)(cid:21) (A.3b)where the derivatives of R τab can be obtained analytically. Appendix B Lattice summations of the DFTB γ -function The phase-modulated lattice sum of the γ function is given by: e γ IJ ( q ) ≡ X R e i qR γ IJ ( R + u J − u I ) (B.1)where u X is the position vector of atom X , as defined in the original unit cell.Adding and subtracting the Coulomb potential, the expression for γ -function can be written as: γ IJ ( r ) = (cid:18) γ IJ ( r ) − r (cid:19) + 1 r (B.2)Since γ IJ ( r ) → /r as r → ∞ , the first term in this expression is short-ranged, making the lattice summationover it straightforward. This is not the case for the second term, so here we use the Ewald summationtechnique, i.e., we split it into a short-ranged and a long-ranged part:1 r = erfc( αr ) r + erf( αr ) r (B.3)where α is an arbitrary positive real number. The short-ranged term can be added to the first term on theRHS of (B.2), whereas the (generalized) Poisson summation formula can be used for the long-ranged part. Theexpression for e γ ( q ) then becomes: e γ IJ ( q ) = X R e i qR γ αIJ ( R + u J − u I ) + X G e V α ( G + q ; u I − u J ) e i ( G + q )( u I − u J ) (B.4)The second sum here runs over all reciprocal vectors G , while α can be chosen to ensure good convergence ofboth sums. γ αIJ ( r ) is a short-ranged function, given by: γ αIJ ( r ) ≡ γ IJ ( r ) − erf( αr ) r + 2 α √ π δ r =0 (B.5)while the expression for e V α ( k ; r ) is more complicated, as it depends on the dimension of the underlying latticeand on whether k k k is finite or not; see table 3 for details.Lattice summations involving the derivatives of the γ -function ((19) and (20)) can be evaluated in a similarmanner. 11able 3: e V α ( k ; r ) for different dimensions. Ω is the measure of the underlying unit cell (i.e., volume, areaand length for three-, two- and one-dimensional systems, respectively). For the two-dimensional case, z is thecomponent of r perpendicular to the direction of the periodicity. For the one-dimensional case, ρ ≡ √ x + y ,where x and y are components of r perpendicular to the direction of the periodicity, Γ( u, v ) is the upperincomplete gamma-function and γ E is the Euler-Mascheroni constant.dim e V α ( k = ; r ) e V α ( k = ; r )3 4 π Ω e − k / α k π Ω 1 k (cid:20) e − kz erfc (cid:18) k α − αz (cid:19) + e kz erfc (cid:18) k α + αz (cid:19)(cid:21) π Ω (cid:20) z erf( αz ) + e − α z α √ π (cid:21) (cid:20) δ ρ =0 ∞ X n =0 ( − n n n ! ( kρ ) n Γ (cid:16) − n, k α (cid:17) + δ ρ =0 Γ (cid:16) , k α (cid:17)(cid:21) − h γ E + Γ(0 , α ρ ) + log( α ρ ) i References [1] Porezag, D.; Frauenheim, T.; K¨ohler, T.; Seifert, G.; Kaschner, R.
Phys. Rev. B , , 12947–12957.[2] Seifert, G.; Porezag, D.; Frauenheim, T. International Journal of Quantum Chemistry , , 185–192.[3] Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B , , 7260–7268.[4] Oliveira, A. F.; Seifert, G.; Heine, T.; Duarte, H. A. J. Braz. Chem. Soc. , , 1193–1205.[5] Gaus, M.; Goez, A.; Elstner, M. Journal of Chemical Theory and Computation , , 338–354.[6] Witek, H. A.; Irle, S.; Morokuma, K. Journal of Chemical Physics , , 5163–5170.[7] Nishimoto, Y.; Irle, S. Chemical Physics Letters , , 317–321.[8] Giannozzi, P.; de Gironcoli, S.; Pavone, P.; Baroni, S. Physical Review B , , 7231–7242.[9] Savrasov, S. Y. Physical Review Letters , , 2819–2822.[10] SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, AMS DFTB 2018. 2018; .[11] Koskinen, P.; M¨akinen, V. Computational Materials Science , , 237–253.[12] Mulliken, R. The Journal of Chemical Physics , , 2343–2346.[13] Visvesvara, R.; Rao, R. Signals and Systems ; Prentice-Hall Of India Pvt. Limited, 2008.[14] Gonze, X.
Physical Review A , , 1096–1114.[15] Ewald, P. P. Annalen der Physik , , 253–287.[16] Sternheimer, R. M. Phys. Rev. , , 951–968.[17] Shcherbakov, M. V.; Brebels, A.; Shcherbakova, N. L.; Tyukov, A. P.; Janovsky, T. A.; Kamaev, V. A. World Applied Sciences Journal , , 171–176.[18] Maradudin, A.; Voske, S. Reviews of Modern Physics , .[19] Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; Van Lenthe, E.; Witek, H. A.; Heine, T. Journal of Chemical Theory and Computation , , 4006–4017.[20] Oliveira, A. F.; Philipsen, P.; Heine, T. Journal of Chemical Theory and Computation , 5209–5218.[21] Rauls, E.; Elsner, J.; Gutierrez, R.; Frauenheim, T.
Solid State Communications , , 459–464.1222] Moreira, N. H.; Dolgonos, G.; Aradi, B.; da Rosa, A. L.; Frauenheim, T. Journal of Chemical Theory andComputation , , 605–614, PMID: 26610226.[23] Markussen, T.; Palsgaard, M.; Stradi, D.; Gunst, T.; Brandbyge, M.; Stokbro, K. Phys. Rev. B , ,245210.[24] Slater, J. C.; Koster, G. F. Physical Review , , 1498–1524.[25] Martial, M. Journal of Physics A: Mathematical and Theoretical , , 425002.[26] Porto, M. Journal of Physics A: Mathematical and General ,33