Wavelets for Density-Functional Theory and Post-Density-Functional-Theory Calculations
Bhaarathi Natarajan, Mark E. Casida, Luigi Genovese, Thierry Deutsch
aa r X i v : . [ c ond - m a t . o t h e r] O c t In: Theoretical and Computational Methods in ModernDensity Functional TheoryEditor: A.K. Roy, pp. 1-45 ISBN 0000000000c (cid:13)
Chapter 1 W AVE L ETS FOR D E NSI TY -F UNCTIONAL T HE ORYAND P OST -D ENSITY -F UNCTIONAL -T HEORY C AL CULATIONS
Bhaarathi Natarajan a,b ∗ Mark E. Casida a † Luigi Genovese b ‡ Thierry Deutsch b § a Laboratoire de Chimie Th´eorique,D´epartement de Chimie Mol´ecularie (DCM, UMR CNRS/UJF 5250),Institut de Chimie Mol´eculaire de Grenoble (ICMG, FR2607),Universit´e Joseph Fourier (Grenoble I),301 rue de la Chimie, BP 53,F-38041 Grenoble Cedex 9, France b UMR-E CEA/UJF-Grenoble 1,INAC, Grenoble, F-38054, FranceOctober 24, 2011 ∗ E-mail address: [email protected] † E-mail address: [email protected] ‡ E-mail address: [email protected] § E-mail address: [email protected]
Natarajan, Genovese, Casida, and Deutsch
PACS
Keywords:
Wavelets, density-functional theory, time-dependent density-functional the-ory, linear-response time-dependent density-functional theory, orbital energies, electronicexcitation energies.
Abstract
We give a fairly comprehensive review of wavelets and of their application todensity-functional theory (DFT) and to our recent application of a wavelet-based ver-sion of linear-response time-dependent DFT (LR-TD-DFT). Our intended audience isquantum chemists and theoretical solid-state and chemical physicists. Wavelets are aFourier-transform-like approach which developed primarily in the latter half of the lastcentury and which was rapidly adapted by engineers in the 1990s because of its advan-tages compared to standard Fourier transform techniques for multiresolution problemswith complicated boundary conditions. High performance computing wavelet codesnow also exist for DFT applications in quantum chemistry and solid-state physics, no-tably the B IG DFT code described in this chapter. After briefly describing the basicequations of DFT and LR-TD-DFT, we discuss how they are solved in B IG DFT andpresent new results on the small test molecule carbon monoxide to show how B IG DFTresults compare against those obtained with the quantum chemistry gaussian-type or-bital (GTO) based code DE M ON K . In general, the two programs give essentially thesame orbital energies, but the wavelet basis of B IG DFT converges to the basis set limitmuch more rapidly than does the GTO basis set of DE M ON K . Wavelet-based LR-TD-DFT is still in its infancy, but our calculations confirm the feasibility of implementingLR-TD-DFT in a wavelet-based code. Contents
1. Introduction 32. Wavelet Theory 7
3. Density Functional Theory 164. Time-Dependent Density Functional Theory 175. Krylov Space Methods 196. Numerical Implementation of DFT in B IG DFT B IG DFT and TD-DFT 28
8. Results 30 DE M ON K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 308.1.2. B IG DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 318.2. Orbital Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 328.3. Excitation Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 338.4. Oscillator Strengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
9. Conclusion 37A List of Abbreviations 38
1. Introduction
The broad meaning of “adaptivity” is the capacity to make something work better by al-ternation, modification, or remodeling. Concepts of adaptivity have found widespread usein quantum chemistry, ranging from the construction of Gaussian-type orbital (GTO) ba-sis sets, see e.g., the development of correlation consistent bases [1, 2, 3], to linear scalingmethods in density functional theory (DFT) [4, 5, 6, 7, 8, 9], selective configuration interac-tion (CI) methods [10, 11] and local correlation methods based on many-body perturbationtheory or coupled cluster (CC) theory [12, 13]. This chapter is about a specific adaptive tool,namely wavelets as an adaptive basis set for DFT calculations which can be automaticallyplaced when and where needed to handle multiresolution problems with difficult boundaryconditions.Let us take a moment to contrast the wavelet concept of adaptivity with other types ofadaptivity. In other contexts, the adaptive procedure is typically based on a combination ofphysical insights together with empirical evidence from numerical simulations. A rigorousmathematical justification is usually missing. This may not be surprising: Familiar con-cepts lose a lot of their original power if one tries to put them in a rigorous mathematicalframework. Therefore, we will not shoulder the monumental and perhaps questionable taskof providing a rigorous mathematical analysis of all the adaptive approaches used nowadaysthroughout quantum chemistry. Instead we will concentrate on the mathematical analysisof a particular electronic structure method which lends itself to a rigorous mathematicalanalysis and application of adaptivity. In contrast with other adaptive methods, multires-olution analysis (MRA) with wavelets can be regarded as an additive subspace correctionand their wavelet representations have a naturally built-in adaptivity which comes through Natarajan, Genovese, Casida, and Deutschtheir ability to express directly and separate components of the desirable functions livingon different scales.This combined with the fact that many operators and their inverses have nearly sparserepresentations in wavelet coordinates may eventually lead to very efficient schemes thatrely on the following principle: Keep the computational work proportional to the numberof significant coefficients in the wavelet expansions of the searched solution. As there area lot of different wavelet bases with different properties (length of support, number of van-ishing moments, symmetry, etc.) in each concrete case we can choose the basis that is mostappropriate for the intrinsic complexity of the sought-after solution. This fact makes thewavelet-based schemes a very sophisticated and powerful tool for compact representationsof rather complicated functions. The expected success of wavelet transforms for solvingelectronic structure problems in quantum mechanics are due to three important properties:(a) the ability to choose a basis set providing good resolution where it is needed, in thosecases where the potential energy varies rapidly in some regions of space, and less in others;(b) economical matrix calculations due to their sparse and banded nature; and (c) the abil-ity to use orthonormal wavelets, thus simplifying the eigenvalue problem. Of course, thismight lead to adaptive methods which are fully competitive from a practical point of view,for example, working with a systematic basis instead of GTO bases requires from the on-set larger basis sets and the benefit of systematic improvement might be a distant prospect.However, we have the more realistic prospect that our rigorous analysis provides new andhopefully enlightening perspectives on standard adaptive methods, which we reckon cannotbe obtained in another way.On the otherhand advances in computational technology opened up new opportunitiesin quantum mechanical calculation of various electronic structures, like molecules, crystals,surfaces, mesoscopic systems, etc. The calculations can only be carried out either for verylimited systems or with restricted models, because of their great demand of computationaland data storage resources. Independent particle approximations, like the Hartree-Fockbased [14, 15, 16, 17] algorithms with single determinant wave functions, leave out theelectron correlation and need operation and storage capacity of order N , if N is the totalnumber of electrons in the system. If inclusion of the electron correlation is necessary, CI orCC methods can be applied, with very high demand of computational resources ( O ( N ) to O ( N !) ). An alternative way is to use MBPT. The second order perturbation calculations canbe carried out within quite reasonable time and resource limits, but the results are usuallyunsatisfactory, they just show the tendencies, while the 4th order MBPT needs O ( N ) to O ( N ) operations. All these algorithms use the N -electron wave function as a basicquantity.Another branch of methods use electron density as the primary entity. Pioneers of thistrend, like Thomas [18], Fermi [19, 20], Frenkel [21] and Sommerfeld [22] developed thestatistical theory of atoms and the local density approximation (LDA). The space aroundthe nuclei is separated into small regions, where the atomic potential is approximated asa constant, and the electrons are modeled as a free electron gas of Fermi-Dirac statistics[23, 24, 21]. Dirac included electron correlation [25], which improved the results. Af-ter the Hohenberg–Kohn theorems had appeared [26], and Kohn and Sham had offered apractically applicable method [27] based on their work, many scientists were motivated towork on the theory, and DFT developed into one of the most powerful electronic structureavelets for DFT and Post-DFT Calculations 5methods.Despite the success of density functional theory, it has some drawbacks. The exactformula of the exchange-correlation potential is not known, thus chemical intuition andmeasured data are necessary in order to approximate it, and the kinetic energy functional ishard to calculate. Powerful approximating formulas are available (see, e.g. [28]), like theThomas–Fermi functional based gradient and generalized gradient expansions, where theenergy functionals are expressed as a power series of the gradient of the density (the firstsuch suggestion was [29].)Considering the historical development of sophisticated N -electron methods, a typicaltrend can be observed. Starting with a very simple model, new details are introduced inorder to improve the results. This scheme is followed in the linear muffin tin orbital method(LMTO) [30] where the interatomic regions is replaced by the spherical orbital of an atomicpotential around the nuclei. Similarly, the linearized augmented plane wave method (APW)[31] and the plane wave pseudopotential approach [32] describe the details of the crystalpotential differently in different spatial regions. Although they are rather successful, forapplying any of these models, chemical intuition is needed, free parameters, like the radiusof the bordering sphere between the two types of potentials, and the boundary conditionshave to be set. A systematic method, which can handle the different behaviors of the elec-tron structures at different spatial domains, or either at different length scale [33], is thelongterm requirement of any physical chemists.Multiresolution or wavelet analysis, this rapidly developing branch of the applied math-ematics, is exactly the tool for statisfy all the need of any chemical physicists/physcialchemists. From mathematical point of view wavelet analysis is a theory of a special kindof Hilbert space basis sets. Basis sets are commonly used in all electron structure calcu-lations, as the wave function is usually expanded as linear combination of some kind ofbasis functions. Thus the operator eigenvalue problem is reduced to an algebraic matrixeigenvector problem. The resulting algebraic equations are easier to solve, well knownalgorithms and subroutine libraries are available, however, the difficulty of choosing theproper basis set arises. If linear combination of atomic orbitals (LCAO) is used, the atomicbasis functions are Slater or Gaussian-type of functions [34, 35], the selection of atomicorbitals needs chemical intuition, which is a result of long time’s experience, and can not bealgorithmized. Both basis sets are non-orthogonal, and lack the explicit convergence prop-erties [36]. Moreover, calculation of operator matrix elements with Slater-type orbitals iscomplicated, their integrals have to be treated numerically. Although integrals of Gaussianfunctions are analytically known, the Gaussisn-type basis does not reflect the nuclear cuspcondition of Kato [37], which reflects on singularities of the N -electron wave function inthe presence of Coulomb-like potentials. Since then it turned out that for high precision nu-merical calculations it is essential to satisfy these requirements. However, while the nuclearcusp condition is relatively easy to fulfill by Slater-type orbitals (STO), the electron-electroncusp is extremely hard to represent. In general, GTO-based/STO-based DFT codes givesreliable results with a relatively small number of basis functions, making them optimal forlarge scale computations where high accuracy is less crucial. On the other hand there is noconsistent way to extend these basis sets and thereby converge the results with respect tothe size of the basis.The second type of basis set covers the system-independent functions such as plane Natarajan, Genovese, Casida, and Deutschwaves [32] or wavelets [38]. The main advantable of these basis sets, is that their size canbe systematically increased until the result of the calculation has converged, and are gen-erally considered to be more accurate than the former type. The number of basis functionrequired to obtain convergence is normally so large that direct solution of the matrix eigen-value problem within the entire basis space is not possible. Instead one has to use iterativemethods to determine the lowest (occpied) part of the spectrum [32]. In solid state physics,where more or less periodic systems are studied, choosing plane wave basis sets is ratherusual. These basis functions are system independent and easily computable, but the resultsare not always convincing and the number of necessary basis functions is almost untreat-able. (Theoretically, plane waves could also be used for describing molecules, since thetwo-electron integrals and the expectation values are connected to the Fourier transform,thus they are easily computable, and this could balance the large number of necessary basisfunctions.) The reason, why so many plane waves are needed is that the wave functionsaround the nuclei need very high frequency terms, i.e. high resolution level, for repro-ducing the nuclear cusps. In the framework of Fourier analysis, the whole space has tobe expanded at the same resolution, despite that in most of the space low frequency termswould be sufficient.Fully-numerical “basis-set free” Hartree-Fock (HF) calculations of atoms have beenknown since the 1960s (Vol. 1, pp. 322-326 and Vol. 2, pp. 15-30 of Ref. [39] andRefs. [40, 41, 42, 43]) and have proven helpful in constructing efficient finite basis setsfor molecular calculations. In the late 1980s, Axel Becke used a fully-numerical density-functional theory (DFT) program for diatomics to show that many of the problems of DFTcalculations at that time were due not to the functionals used, but rather numerical artifactsof the DFT programs of the 1970s [44].) Since that time, fully-numerical DFT codes havebeen implemented for polyatomic molecules using the finite element method (FEM), withP ARSEC from the chemists point of view or OCTOPUS from the view of physicsts beinga notable example.B IG DFT the pseudo potential code for bigger systems based as it is on traditionalHohenberg-Kohn-Sham DFT [26, 27], could only calculate ground-state properties withan eye to order-N DFT. As a step to increase the feasibility of the code we formulatedthe wavelet-based linear-response time-dependent density-functional theory (TD-DFT) andhere we support our first implementation for calculating electronic excitation spectra [45].Electronic excitation spectra can be calculated from TD-DFT [46] using time-dependentlinear response (LR) theory [47, 48]. Casida formulated LR-TD-DFT (often just referedto as TD-DFT) so as to resemble the linear-response time-dependent HF equations alreadyfamiliar to quantum chemists [48]. That method was then rapidly implemented in a largenumber of electronic structure codes in quantum chemistry, beginning with the DE M ON family of programs [49] and the T URBO M OL program [50]. Among the programs that im-plemented “Casida’s equations” early on was the FEM DFT program PARSEC [51] andalso be found in the FEM DFT program O CTOPUS [52]. See Ref. [53] for a recent FEMimplementation of TD-DFT. Since a wavelet-based program offers certain advantages overthese other FEM DFT programs, it was deemed important to also implement LR-TD-DFTin B IG DFT.In the next section we give a detailed description of the idea behind the multiresolutionanalysis and wavelets, with a historical note. Sec. 3. and Sec. 4., briefly presenting the the-avelets for DFT and Post-DFT Calculations 7oretical introduction to DFT and TD-DFT, and Sec. 5., talks about the well-known Krylovspace methods for solving eigenvalue equations involved in our implementation. Sec. 6.and Sec. 7., gives the numerical implementation of DFT and how we have implementedTD-DFT from the aspects of theoretical and algorithmic point of view on wavelets basedpseudopotential electronic structure code B IG DFT, and in Section 8. we give the resultsof detailed comparisons between TD-DFT excitation spectra calculated with B IG DFT andwith the implementation of Casida’s equations in the GTO-based program DE M ON K . Theconclusion were drawn for future applications in the field of chemistry and some of theother problems are reviewed to draw chemists’ greater attention to wavelets and to gainmore benifits from using wavelet technique.
2. Wavelet Theory
The mathematics of wavelets is a fairly new technique, it can generally be used whereone traditionally uses Fourier techniques. They incorporate the feature of having multiplescales, so very different resolutions can be used in different parts of space in a mathemat-ically rigorous manner. This matches many systems in nature well, for example moleculewhere the atomic orbitals are very detailed close to the cores, while they only vary slowlybetween them. Wavelet analysis can quite generally be viewed as a local Fourier analysis.From the wavelet expansion, or wavelet spectrum, of a function, f , it can be inferred notonly how fast f varies, i.e. which frequencies it contains, but also where in space a givenfrequency is located. This property has important applications in both data compression,signal/image processing and noise reduction [54]. Wavelet methods are also employed forsolving partial differential equations [55, 56], and in relation to electronic structure methodsa complete DFT program based on interpolating wavelets has been developed [57]. Most historical versions of wavelet theory however, despite their source’s perspective, beginwith Joseph Fourier. In 1807, a French mathematician, Joseph Fourier, discovered that allperiodic functions could be expressed as a weighted sum of basic trigonometric functions.His ideas faced much criticism from Lagrange, Legendre and Laplace for lack of mathe-matical rigor and generality, and his papers were denied publication. It took Fourier over15 years to convince them and publish his results. Over the next 150 years his ideas wereexpanded and generalized for non-periodic functions and discrete time sequences. The fastFourier transform algorithm, devised by Cooley and Tukey in 1965 placed the crown onFourier transform, making it the king of all transforms. Since then Fourier transforms havebeen the most widely used, and often misused, mathematical tool in not only electrical en-gineering, but in many disciplines requiring function analysis. This crown however, wasabout to change hands. Following a remarkably similar history of development, the wavelettransform is rapidly gaining popularity and recognition.The first mention of wavelets was in a 1909 dissertation by Hungarian mathematicianAlfred Haar. Haar’s work was not necessarily about wavelets, as “wavelets” would notappear in their current form until the late 1980s. Specifically, Haar focused on orthogonalfunction systems, and proposed an orthogonal basis, now known as the Haar wavelet basis, Natarajan, Genovese, Casida, and Deutschin which functions were to be transformed by two basis functions. One basis function isconstant on a fixed interval, and is known as the scaling function. The other basis functionis a step function that contains exactly one zero–crossing (vanishing moment) over a fixedinterval (more on this later).The next major contribution to wavelet theory was from a 1930s French scientist PaulPierre L´evy. More correctly, L´evy’s contribution was less of a contribution and more ofa validation. While studying the ins and outs of Brownian motion in the years followingHaar’s publication, L´evy discovered that a scale–varied Haar basis produced a more accu-rate representation of Brownian motion than did the Fourier basis. L´evy, being more of aphysicist than mathematician, moved on to make large contributions to our understandingof stochastic processes.Contributions to wavelet theory between the 1930s and 1970s were slight. Most impor-tantly, the windowed Fourier transform was developed, with the largest contribution beingmade by another Hungarian named Dennis Gabor. The next major advancement in wavelettheory is considered to be that of Jean Morlet in the late 1970s.Morlet, a French geophysicist working with windowed Fourier transforms, discoveredthat fixing frequency and stretching or compressing (scaling) the time window was a moreuseful approach than varying frequency and fixing scale. Furthermore, these windows wereall generated by dilation or compression of a prototype Gaussian. These window functionshad compact support both in time and in frequency (since the Fourier transform of a Gaus-sian is also a Gaussian.) Due to the small and oscillatory nature of these window functions,Morlet named his functions as “wavelets of constant shape”. In 1981, Morlet worked withCroatian–French physicist Alex Grossman on the idea that a function could be transformedby a wavelet basis and transformed back without loss of information, thereby outlining thewavelet transformation. It is of note that Morlet initially developed his ideas with nothingmore than a handheld calculator.In 1986, St´ephane Mallat noticed a publication by Yves Meyer that built on the conceptsof Morlet and Grossman. Mallat sought Meyer’s consult, and the result of said consult wasMallat’s publication of multiresolution analysis. Mallat’s MRA connected wavelet trans-formations with the field of digital signal processing. Specifically, Mallat developed thewavelet transformation as a multiresolution approximation produced by a pair of digital fil-ters. The scaling and wavelet functions that constitute a wavelet basis are represented bya pair of finite impulse response filters, and the wavelet transformation is computed as theconvolution of these filters with the input function. The importance of Mallat’s contribu-tion cannot be overstated. Without the fast computational means of wavelet transformationprovided by the MRA, wavelets, undoubtedly, would not be the effective and widely usedsignal processing tools that they are today.In 1988, a student of Alex Grossman, named Ingrid Daubechies, combined the ideas ofMorlet, Grossman, Mallat, and Meyer by developing the first family of wavelets as they areknown today. Named the Daubechies wavelets, the family consists of 8 separate wavelet andscaling functions (more on this later). With the development of pair Daubechies wavelet andscaling functions is orthogonal, continuous, regular, and compactly supported, the founda-tions of the modern wavelet theory were laid. The last ten years mostly witnessed a searchfor other wavelets with different properties and modifications of the MRA algorithms. In1992, Albert Cohen, Jean Feauveau and Daubechies constructed the compactly supportedavelets for DFT and Post-DFT Calculations 9biorthogonal wavelets, which are preferred by many researchers over the orthonormal basisfunctions, whereas R. Coifman, Meyer and Victor Wickerhauser developed wavelet pack-ers, a natural extension of MRA.
A suitable gateway to the theory of wavelets is through the idea of MRA. A detailed de-scription of MRAs can be found in Keinert [58], from which a brief summary of the keyissues are given in the following.A multiresolution analysis is an infinite nested sequence of subspaces L ( R ) V j ⊂ V j ⊂ ... ⊂ V nj ⊂ ... (1)with the following properties • V ∞ j is dense in L • f ( x ) ∈ V nj ⇐⇒ f (2 x ) ∈ V n +1 j ≤ n ≤ ∞• f ( x ) ∈ V nj ⇐⇒ f ( x − − n l ) ∈ V nj ≤ l ≤ (2 n − • There exists a function vector ϕ of length j + 1 in L such that { ϕ j ( x ) : ≤ k ≤ j } forms a basis for V j . This means that if we can construct a basis of V j , which consists of only j +1 functions,we can construct a basis of any space V nj , by simple compression (by a factor of n ), andtranslations (to all grid points at scale n ), of the original j + 1 functions, and by increasingthe scale n , we are approaching a complete basis of L . Since V nj ⊂ V n +1 j the basisfunctions of V nj can be expanded in the basis of V n +1 j ϕ nl ( x ) def = 2 n/ ϕ (2 n x − l ) = X l h ( l ) ϕ n +1 l ( x ) . (2)where h ( l ) s are the so-called filter matrix that describes the transformation between differentspaces V nj .The MRA is called orthogonal if h ϕ n ( x ) , ϕ nl ( x ) i = δ l I j +1 , (3)where I j +1 is the ( j + 1) × ( j + 1) unit matrix, and j + 1 is the length of the functionvector. The orthogonality condition means that the functions are orthogonal both withinone function vector and through all possible translations on one scale, but not through thedifferent scales.Complementary to the nested sequence of subspaces V nj , we can define another seriesof spaces W nj that complements V nj in V n +1 j V n +1 j = V nj ⊕ W nj (4)0 Natarajan, Genovese, Casida, and Deutschwhere there exists another function vector φ of length j + 1 that, with all its translations onscale n forms a basis for W nj . Analogously to Eq. (2) the function vector can be expandedin the basis of V n +1 j φ nl ( x ) def = 2 n/ φ (2 n x − l ) = X l g ( l ) φ n +1 l ( x ) . (5)with filter matrices g ( l ) . The functions φ also fulfill the same orthogonality condition asEq. (3), and if we combine Eq. (1) and Eq. (4) we see that they must be orthogonal withrespect to different scales. Using Eq. (4) recursively we obtain V nj = V j ⊕ W j ⊕ W j ⊕ ... ⊕ W n − j . (6)which will prove to be an important relation. There are many ways to choose the basis functions ϕ and φ (which define the spannedspaces V nj and W nj ), and there have been constructed functions with a variety of properties,and we should choose the wavelet family that best suits the needs of the problem we aretrying to solve. (Wavelets are often denoted by ψ in the literature but the choice has beenmade here to denote them by φ so as to avoid confusion with the Kohn-Sham orbitals.)Otherwise, we could start from scratch and construct the new family, one that is custum-made for the problem at hand. Of course, this is not a trivial task, and it might prove moreefficient to use an existing family, even though its properties are not right on cue.There is a one-to-one correspondence between the basis functions ϕ and φ , and the filtermatrices h ( l ) and g ( l ) used in the two-scale relation equations Eq. (2) and Eq. (5), and mostwell-known wavelet families are defined only through their filter coefficients.In the following we are taking a different, more intuitive approach, for defining the scaling space V nj as the space of piecewise polynomial functions V nj def = { f : all polynomials of degree ≤ j on the interval (2 − n l, − n ( l + 1)) for ≤ l < n , f vanishes elsewhere } . (7)It is quite obvious that one polynomial of degree j on the interval [0 , can be exactlyreproduced by two polynomials of degree j , one on the interval [0 , ] and the other onthe interval [ , . The spaces V nj hence fulfills the MRA condition Eq. (1), and if thepolynomial basis is chosen to be orthogonal, the V nj constitutes an orthogonal basis. The basic wavelet ideas that we need can be easily explained using Haar wavelets [59].These are simply the box functions shown in Fig. 2. We begin with a compact “motheravelets for DFT and Post-DFT Calculations 11Figure 1. Wavelets (bottom) and scaling function (top).scaling function,” in this case the Haar function, ϕ ( x ) = x >
11 ; 0 < x <
10 ; x < . (8)Translations, { ϕ i ( x ) = ϕ ( x − i ) } , of this mother function produces a crude basis set. Itsrelation to the grid of integers is obvious. Successively more refined basis sets may begenerated by repeated application of the scaling operation consisting of contracting thefunctions to half their size in the x direction. The k th generation of scaling function isgiven by n ϕ ( k ) i ( x ) = ϕ (2 k x − i ) o . Each generation has a fixed resolution related to anunderlying grid with the same resolution. Let us now try to construct a multiresolution basisset. This is accomplished by (say) beginning with the third generation wavelets and takingsums and differences of adjacent functions until the eight third generation scaling functionshave been replaced with the eight wavelets shown in Fig. 3. Notice how each generation ofdaughter wavelets is related to the mother wavelet by scaling, φ ( k ) i ( x ) = φ (2 k x − i ) . Noticealso how the mother and two generations of daughter wavelets plus the mother scalingfunction (occasionally refered to as the “father wavelet”) constitute a multiresolution basisset equivalent to the original third generation scaling basis set. Thus an arbitrary function, f ( x ) , expressible in the original scaling basis, f ( x ) = X i =1 ϕ (3) i ( x ) s (3) i , (9)has the wavelet transform, f ( x ) = ϕ (0)0 ( x ) s (0)0 + φ (0)0 ( x ) d (0)0 + X i =0 , φ (1) i ( x ) d (1) i + X i =0 , φ (2) i ( x ) d (2) i . (10)Since the basis set is multiresolution, we may choose to add more grid points in some regionof space and go locally to higher order wavelet expansions. It is also not always necessaryto carry out a full wavelet transform, but rather it may be useful to just carry out a partialtransform giving a linear combination of wavelets with several scaling functions at a time.2 Natarajan, Genovese, Casida, and DeutschFigure 2. Haar scaling functions.Figure 3. Haar scaling functions and the corresponding wavelets.avelets for DFT and Post-DFT Calculations 13The extension to three dimensions is accomplished by using products of one-dimensionalscaling functions and wavelets. Haar wavelets are just one type of wavelet basis set. Ithappens to be pedagogically useful but is not particularly useful for computations. The wavelet space W nj is defined, according to Eq. (4), as the orthogonal complement of V nj in V n +1 j . The wavelet basis functions of W nj are hence piece-wise polynomials of degree ≤ j on each of the two intervals on scale n + 1 that overlaps with one interval on scale n .These piece-wise polynomials are then made orthogonal to a basis of V nj and to each other.The construction of the wavelet basis follows exactly [60] where a simple Gram-Schmidtorthogonalization were employed to construct a basis that met the necessary orthogonalityconditions.One important property of the wavelet basis is the number of vanishing moments. The j -th continuous moment of a function φ is defined as the integrals µ j def = Z x j φ ( x ) dx , (11)and the function φ has M vanishing moments if µ j = 0 , k = 0 , ..., M − The vanishing momenets of the wavelet functions gives information on the approxima-tion order of the scaling functions. If the wavelet function φ has M vanishing moments,any polynomial of order ≤ M − can be exactly reproduced by the scaling function ϕ ,and the error in representing an arbitrary function in the scaling basis is of M -th order. Byconstruction, x i is in the space V j for ≤ i ≤ j , and since W j ⊥ V j , the first k + 1 moments of φ j must vanish. The construction of the scaling functions is quite straightforward, j +1 suitable polynomialsare chosen to span any polynomial of degree j on the unit interval. The total basis for V nj is then obtained by appropriate dilation and translation of these functions. Of course,any polynomial basis can be used, the simplest of them the standard basis { , x, ..., x j } .However, this basis is not orthogonal on the unit interval. In the following, two choices oforthogonal scaling functions will be presented, and even though they span exactly the samespaces V nj there are some important numerical differences between the two.In order to construct a set of orthogonal polynomials we could proceed in the same man-ner as for the wavelet functions and do a Gram-Schmidt orthogonalization of the standardbasis { , x, ...x j } . If this is done on the interval x ∈ [ − , we end up with the Legendrepolynomials { L k } jk =0 . These functions are usually normalized such that L k (1) = 1 for all j . To make the Legendre scaling functions ϕ Lk we transform the Legendre polynomials tothe interval x ∈ [0 , , and L normalize ϕ Lk ( x ) = √ k + 1 L k (2 x − , x ∈ [0 , . (12)4 Natarajan, Genovese, Casida, and DeutschThe basis for the space V nj is then made by proper dilation and translation of ϕ Lk . Alpert etal. [60] presented an alternative set of scaling functions with interpolating properties. These interpolating scaling functions ϕ Ik are based on the Legendre scaling functions { ϕ Lk } jk =0 ,and the roots { y k } jk =0 and weights { w k } jk =0 of the Gauss-Legendre quadrature of order j + 1 , and are organized in the linear combinations ϕ Ik ( x ) = √ w k j p X i =0 ϕ Li ( y k ) ϕ Li ( x ) , x ∈ [0 , , (13)Again the basis of V nj is made by dilation and translation of φ Ik . The construction of ϕ Ik gives them the interpolating property ϕ Ik ( y i ) = δ ki √ w i . (14)which will prove important for numerical efficiency.A detailed discussion on the properties of interpolating wavelets can be found inDonoho [61], but the case of interpolating wavelets is somewhat different. An importantproperty of interpolating wavelets is the smoothness of any function represented in thisbasis. This property stems from general Lagrange interpolation. In the wavelet case theinterpolating property applies within one scaling function vector only, which means thatfunctions represented in this basis can be discontinous in any merging point between thedifferent translations on any scale. Since the general introduction to wavelets were already made, we will now concencentrateour description on the level 3 interpolating scaling function (ISF) introduced by Deslauriersand Dubuc, and described in detail in Ref. [62]. Its main advantage is that it is fast andeasy to perform nonlinear operations on functions represented in this basis, as long as theoperation is local in shape. It also represents 3rd order polynomials exactly which meansthat it behaves very smoothly.We introduced the projection operator P n that projects an arbitrary function f ( x ) ontothe basis { ϕ nj,l } of the scaling space V n (in the remaining of this text the subscript k ofthe scaling and wavelet spaces will be omitted, and it will always be assumed that we aredealing with a k th order polynomial basis.) f ( x ) ≈ P n f ( x ) def = f n ( x ) = n − X l =0 k X j =0 s n,fj,l ϕ nj,l ( x ) , (15)where the expansion coefficients s n,fj,l , the so-called scaling coefficients, are obtained bythe usual integral s n,fj,l def = h f, ϕ nj,l i = Z f ( x ) ϕ nj,l ( x ) dx , (16)If this approximation turns out to be too crude, we double our basis set by increasing thescale and perform the projection P n +1 . This can be continued until we reach a scale N where we are satisfied with the overall accuracy of f N relative to the true function f .avelets for DFT and Post-DFT Calculations 15In a perfect world, the projection in Eq. (16) could be done exactly, and the accuracy ofthe projection would be independent of the choice of polynomial basis. In the real world theprojections are done with Gauss-Legendre quadrature and the expansion coefficients s n,fj,l of f ( x ) are obtained as s n,fj,l = Z − n ( l +1) − nl f ( x ) ϕ nj,l ( x ) dx = 2 − n/ Z f (2 − n ( x + l )) ϕ j, ( x ) dx ≈ − n/ k q − X q =0 w q f (2 − n ( y q + l )) ϕ j, ( y q ) (17)where { w q } k q − q =0 are the weights and { y q } k q − q =0 the roots of the Legendre polynomial L k q used in k q th order quadrature.By approximating this integral by quadrature we will of course not obtain the exactexpansion coefficients. However, it would be nice if we could obtain the exact coefficientswhenever our basis is flexible enough to reproduce the function exactly, that is if f ( x ) isa polynomial of degree ≤ k . The Legendre quadrature holds a (2 k − rule which statesthat the k -order quadrature is exact whenever the integrand is a polynomial of order k − .By choosing k q = k + 1 order quadrature we will obtain the exact coefficient whenever f ( x ) is a polynomail of degree ≤ ( k + 1) when projecting on the basis of k -order Legendrepolynomials.In the multidimensional case the expansion coefficients are given by multidimensionalquadrature s nfjl = 2 − dn/ k X q =0 k X q =0 ... k X q d =0 f (2 − n ( y q + l ))Π di =1 w q i ϕ j p , ( y q i ) , (18)using the following notation for the vector of quadrature roots y q def = ( y q , y q , ..., y q d ) , (19)This quadrature is not very efficient in multiple dimensions since the number of terms scalesas ( k + 1) d . However, if the function f is separable and can be written f ( x , x , ..., x d ) = f ( x ) f ( x ) ...f d ( x d ) , Eq. (18) can be simplified to s nfjl = 2 − dn/ Π di =1 k X q i =0 f i (2 − n ( y q i + l i )) w q i ϕ j i , ( y q i ) , (20)which is a product of small summations and scales only as d ( k + 1) .The Legendre polynomials show very good convergence for polynomial functions f ( x ) ,and are likely to give more accurate projections. However, most interesting functions f ( x ) are not simple polynomials, and the accuracy of the Legendre scaling functions versus ageneral polynomial basis might not be very different.6 Natarajan, Genovese, Casida, and DeutschBy choosing the quadrature order to be k +1 a very important property of the Interpolat-ing scaling functions emerges, stemming from the specific construction of these functionsEq. (13), and the use of the k + 1 order quadrature roots and weights. The interpolatingproperty Eq. (14) inserts a Kronecker delta whenever the scaling function is evaluated in aquadrature root, which is exactly the case in the quadrature sum. This reduces Eq. (17) to s n,fjl = 2 − n/ √ w j f (2 − n ( x j + l )) , (21)which obviously makes the projection k + 1 times more efficient.In multiple dimensions this property becomes even more important, since it effectivelyremoves all the nested summations in Eq. (18) and leaves only one term in the projection s nfjl = f (2 − n ( y j + l ))Π di =1 − n/ √ w j i , (22)This means that in the Interpolating basis the projection is equally effective regardless ofthe separability of the function f .
3. Density Functional Theory
A method to resolve the electronic structure is by using variational principle E [Ψ] = h Ψ | ˆ H | Ψ ih Ψ | Ψ i , (23)Where h Ψ | ˆ H | Ψ i = R d r Ψ ∗ ( r ) ˆ H Ψ( r ) , Ψ denotes the electronic wavefunction and ˆ H theHamiltonian. The energy computed from a guess Ψ is an upper bound to the true groundstate energy E . Full minimization of the functional E [Ψ] will give the true ground state Ψ gs and energy E = E [Ψ gs ] .Density-functional theory states that the many electron problem can be replaced by anequivalent set of self-consistent one-electron equations, the Kohn-Sham equations ˆ hψ σi ( r ) = (cid:18) − ∇ + ˆ v pp ( r ) + ˆ v H ( r ) + ˆ v σxc ( r ) (cid:19) ψ σi ( r ) = ǫ σi ψ σi ( r ) . (24)The eigenfunctions ψ σi are the one-electron wavefunctions that correspond to the minimumof the Kohn-Sham energy functional. In these wavefunctions, i is the orbital index and σ denotes the spin, which can be either up ↑ or down ↓ (spin α or β .)The Hamiltonian ˆ H consists of four different parts: a part related to the kinetic energyof the electrons, the pseudopotentials ˆ v psp , the Hartree potential ˆ v H and the exchange cor-relation potential ˆ v xc . The interaction of the positively charged nuclei with the electrons isdescribed using the pseudopotential ˆ v psp instead of using the full Coulombic potential. Thepseudopotential usually consists of both a local and a non-local part ˆ v psp ( r ) = v loc ( r ) + X l | l i ˆ v l ( r, r ′ ) h l | . (25)avelets for DFT and Post-DFT Calculations 17The Hartree potential ˆ v H describes the interaction between electrons and is given by ˆ v H ( r ) = Z d r ′ ρ ↑ ( r ) + ρ ↓ ( r ′ ) | r − r ′ | . (26)Finally, the exchange correlation potential ˆ v xc describes the nonclassical interaction be-tween the electrons and is given by the functional derivative of an exchange correlationenergy functional ˆ v σxc ( r ) = δE xc ( ρ ↑ , ρ ↓ ) δρ σ ( r ) . (27)In these equations ρ σ is the electron spin density, defined as ρ σ ( r ) = X i n σi | ψ σi ( r ) | , (28)where n σi is the occupation number, i.e. the number of electrons in orbital i . In case of LDA(which we use throughtout this chapter) where there is no longer a distinction between spinup and spin down, orbitals can contain at most two electrons.
4. Time-Dependent Density Functional Theory
This section contains a brief review of the basic formalism of TD-DFT which is alreadywell-known from the literature [46]. The time-dependent single particle Kohn-Sham equa-tions are, (cid:18) − ∇ + v eff [ ρ ]( r , t ) (cid:19) ψ iσ ( r , t ) = i ∂∂t ψ iσ ( r , t ) (29)Here, the wave functions ψ i ( r , t ) and v eff [ ρ ]( r , t ) explicitly depend on time, whereas, v eff [ ρ ]( r , t ) = X a v ion ( r − R a ) + Z ρ [ r ′ , t ] | r − r ′ | d r ′ + v xc [ ρ ]( r , t ) . (30)Using the adiabatic approximation (AA), (which is local in time) v xc [ ρ ]( r , t ) ≈ δE xc [ ρ ] δρ ( r ) δv xc [ ρ ]( r , t ) δρ ( r ′ , t ) ≈ δ ( t − t ′ ) δ E xc [ ρ ] δρ ( r ) δρ ( r ′ ) , (31)and using the LDA, E xc [ ρ ] = Z ρ ( r ) ǫ xc ( ρ ( r )) d r . (32)The method that we wish to use here is Casida’s approach [48]. This section explainsthe deriving equations of linear-response (LR) TD-LDA method.The time-dependent perturbation to the external potential can be written as, δv eff [ ρ ]( r , t ) = δv appl ( r , t ) + δv SCF [ ρ ]( r , t ) (33)8 Natarajan, Genovese, Casida, and Deutschwhere, v SCF [ ρ ]( r , t ) = Z ρ ( r ′ , t ) | r − r ′ | d r ′ + v xc [ ρ ]( r ′ , t ) (34)The LR of the density matrix (DM) can be written in terms of generalised susceptibility χ as below, δ P ijσ ( ω ) = X klτ χ ijσ,klτ ( ω ) δv effklτ ( ω ) (35)After some mathematical steps, one can end-up with the sum-over-states (SOS) representa-tion of χ , χ ijσ,klτ ( ω ) = δ ik δ jl δ στ λ lkτ ω − ( ω lkτ ) (36)where λ lkτ = n lτ − n kτ the difference in occupation numbers and ω lkτ = ǫ kτ − ǫ lτ the difference between the eigenvalues of l th and k th states. In the basis of Kohn-Shamorbitals { ψ iσ } , we can re-write the LR-DM equation as, δ P ijσ ( ω ) = λ jiσ ω − ω jiσ h δv applijσ ( ω ) + δv SCFijσ ( ω ) i (37)Now the term δv SCF is complicated because it itself depends on the response of the DM. δv SCFijσ ( ω ) = X K ijσ,klτ δ P klτ ( ω ) (38)Where, K ijσ,klτ = ∂v SCFijσ ∂ P klτ (39)whose integral form is, K ijσ,klτ = Z Z ψ ∗ iσ ( r ) ψ jσ ( r ) (cid:20) | r − r ′ | + δ E xc [ ρ ] δρ σ ( r ) δρ τ ( r ′ ) (cid:21) ψ kτ ( r ′ ) ψ ∗ lτ ( r ′ ) d r d r ′ (40)If the response is due to a real spin independent external perturbation, δv appl , then we canrestrict ourselves to the real density response and the coupling matrix is symmetric.After some algebra, the real parts of the DM elements ℜ δ P ( ω ) can be given as, X klτ (cid:20) δ ik , δ jl δ στ λ klτ ω klτ ( ω − ω klτ ) − K ijσ,klτ (cid:21) ℜ ( δ P klτ )( ω ) = δv applijσ ( ω ) (41)Here the real part of ℜ δ P σ ( ω ) means the Fourier transform of the real part of ℜ δ P σ ( t ) .Thus the real part of the first-order DM obtained from the solution of the above linear equa-tions gives access to the frequency-dependent polarizabilities. This leads to the followingeigenvalue equation from which the excitation energies and oscillator strengths can be ob-tained. Ω ~F I = ω I ~F I , (42)where, Ω ijσ,klτ = δ ik δ jl δ στ ω klτ + 2 p λ ijσ ω ijσ K ijσ,klτ p λ klτ ω klτ (43)avelets for DFT and Post-DFT Calculations 19Here, the desired excitation energies are equal to ω I and the oscillator strengths f I are ob-tained from the eigenvectors ~F I . The frequency-dependent polarizability is directly realtedto vertical excitation energies, oscillator strength and transition dipole moments µ I , α ( ω ) = X I f I ω I − ω = 23 X I ω I µ I ω I − ω (44)
5. Krylov Space Methods
The methods described in this article involve solving very large eigenvalue problems. Oneof these is the matrix form of the Kohn-Sham orbital equation Eq. (24) while the other is theLR-TD-DFT equation Eq. (42). The first is very large because the wavelet basis set is verylarge while the other is very large because it is of the order of the number of unoccupiedorbitals times the number of unoccupied orbitals on each side. It is important to realize thatspecial methods must be and are being used to solve these very large eigenvalue problems.In particular, B IG DFT make use of the block Davidson variant of the Krylov space methodto solve the Kohn-Sham equation while B IG DFT and most other codes solving the LR-TD-DFT equation Eq. (42) also make use of the the block Davidson method. Krylov methodsand the block Davidson method are briefly described in this section.Krylov space methods may be traced back to a paper in the early 1930s written bythe Russian mathematician Alexei Nikolaevich Krylov. The main idea is that to solve amatrix problem involving a matrix A , it is frequently never actually necessary to constructthe matrix A because iterative solutions only require a reasonable first guess followed byrepeated action of the operator ˆ A on a vector. A number of such methods are known withLanczos diagonalization and the discrete inversion in the iterative subspace (DIIS) [63] asparticlarly well-known examples in theoretical chemical physics. Given a vector ~x , theKrylov space of order r is given by, K r ( A , ~x ) = span (cid:8) ~x, A ~x, A ~x, · · · , A r ~x (cid:9) . (45)The Davidson diagonalization method [64] for solving the matrix eigenvalue problem A ~x = a~x , (46)is deceptively simple. Suppose that we want the lowest eigenvalue and eigenvector and wehave an intial guess vector, ~x (0) . Then we can always write, ~x = ~x (0) + δ~x , (47)is the component of the exact solution which is orthogonal to the intial gues vector. Sim-ple algebra then gives a formula highly reminiscent of Rayleigh-Schr¨odinger perturbationtheory but exact, δ~x = [ Q ( a − A ) Q ] − ( A − a ) ~x (0) , (48)where, Q = − ~x (0) ~x (0) † , (49)0 Natarajan, Genovese, Casida, and Deutschprojects onto the subspace orthogonal to the guess vector. Solving Eq. (48) requires us toovercome two difficulties. The first is that we need a guess for the eigenvalue a , but thisis easily remedied by taking a (0) = ~x (0) † A ~x (0) /~x (0) † ~x (0) and then iterating. The problemgreater problem is to invert the matrix [ Q ( a − A ) Q ] . It actually turns out that a highlyaccurate inversion is not really needed (and actually can cause some problems.) Instead itis better just to replace this difficult inversion with, δ~x = ( a − D ) − ( A − a ) ~x (0) , (50)where D is some diagonal matrix, hence easy to invert. Orthogonalizing δ~x to ~x (0) andnormalizing produces ~x (1) , which is the next basis vector in our Krylov space. Setting upand diagonalizing the 2 × A in this basis and taking the lowest eigenvalue givesus the next estimate a (1) . If application of Eq. (48) is close to zero then we have solvedthe eigenvalue problem Eq. (46), otherwise we have a new δ~x with which to generate ~x (2) and so on and so forth until convergence. The block Davidson method [65] extends theDavidson method to the lowest several eigenvalues and eigenvectors.Davidson diagonalization works well when started from a reasonably good initial guess,otherwise the Lanczos method may be advantageous. One of the most recent incarnations ofthe Lanczos method is the Liouville-Lanczos method for solving the LR-TD-DFT problem[66, 67, 68, 69].
6. Numerical Implementation of DFT in B IG DFT
Computational physics/chemistry is the transformation and implementation of scientifictheory into efficient algorithms which requires both theoretical and experimental skill. Thetransformation of a new theory into an efficient algorithm requires understanding of pro-gramming concepts, mathematical and physical intuition and theoretical insight, whereasthe production of the computer code is much like experimentation, requiring debugging,testing and organisation to yield a highly efficient product. It is also an adaptation of newscientific theory into computer code exploiting the advances in compiler, programming lan-guage and hardware technology. The aim is to afford an algorithm to enable efficient com-putation, portability of code, ease of adaptability and to document the science. To affordsuch an algorithm requires an intuitive understanding of the physics to be implemented,much experimentation with optimisation and debugging of the developing code, a suitablechoice of programming language, as well as a basic overview of the nature of the platformsfor which the code is intended.The Kohn-Sham scheme of DFT greatly reduces the complexity of ground state elec-tronic structure calculations by recasting the many-body problem into a (self-consistent)single-particle problem. For real atomistic systems, however, the KS equations are stilldifficult to solve and further approximate techniques are required. In general it is impor-tant, though, that these approximations can be controlled in such a way that the associatederror does not exceed the error introduced by the xc-functional. While DFT accounts forapproximately 90 % of all quantum chemical calculations being performed, the sometimesunpredictable nature of results and the inability to systematically improve the quality ofcalculation may mean that a place for the conventional correlated techniques remains in theavelets for DFT and Post-DFT Calculations 21quantum chemist’s tool kit. In this work the detailed description of DFT program B IG DFThas been given. B IG DFT [57] has been developed as an European project (FP6-NEST)from 2005 to 2008, and is a wavelet-based pseudopotential implementation of DFT andTD-DFT. For complimentary purpose, Gaussian based quantum chemistry DFT code DE -M ON K [70] is also used but we are not going to discuss the numerical implementation of DE M ON K here and we restrict ourselves to use DE M ON K for revalidating our recent im-plementation of LR-TD-DFT in B IG DFT. However in the following sections, we are goingto recast how the fundamental computational operations were performed in B IG DFT.
Before embarking on our own endeavours, we should make some reference to related work.First, it should be acknowledged that a considerable amount of work has been done alreadyin pursuit of a wavelets in the electronic structure calculations [71, 72, 73, 74, 75, 76].The object of using wavelets as basis set is to associate an expansion coefficient to eachof the piece-wise wavelets. The expansion coefficients are free to vary from one waveletfunction to the next. This feature enables wavelets as highly localised continuous functionsof a fractal nature that have finite supports. The Daubechies wavelets have no availableanalytic forms, and they are not readily available in sampled versions. They are definedeffectively by the associated dilation coefficients. These express a wavelet in high resolutionand a scaling function in the low resolution–which has the same width and which stretchesto zero–as a linear combination of the more densely packed and less dispersed scalingfunctions that form a basis for the two resolution level in combination.The fact that the Daubechies wavelets are known only via their dilation coefficients is noimplediment to the discrete wavelet transform. This transform generates the expansion co-efficients associated with the wavelet decomposition of a data sequence. In this perspective,the dilation coefficients of the wavelets and of the associated scaling functions are nothingbut the coefficients of a pair of quadrature mirror filters that are applied succesively.As described above Daubechies family consists of two fundamental functions: the scal-ing function φ ( x ) and the wavelet ϕ ( x ) (see Fig. 4.) The full basis set can be obtained fromall translations by a certain grid spacing h of the scaling and wavelet functions centered atthe origin. These functions satisfy the fundamental defining (refinement) equations, φ ( x ) = √ m X j =1 − m h j φ (2 x − j ) ,ϕ ( x ) = √ m X j =1 − m g j φ (2 x − j ) . (51)which relate the basis functions on a grid with spacing h and another one with spacing h/ .The coefficients, h j and g j , consitute the so-called “filters” which define the wavelet familyof order m . These coefficients satisfy the relations, P j h j = 1 and g j = ( − j h − j +1 .Eq. (51) is very important since it means that a scaling-function basis defined over a finegrid of spacing h/ may be replaced by combining a scaling-function basis over a coarsegrid of spacing h with a wavelet basis defined over the fine grid of spacing h/ . This thengives us the liberty to begin with a coarse description in terms of scaling functions and2 Natarajan, Genovese, Casida, and Deutschthen add wavelets only where a more refined description is needed. In principle the refinedwavelet description may be further refined by adding higher-order wavelets where needed.However in B IG DFT we restricted ourselves to just two levels: coarse and fine associatedrespectively with scaling functions and wavelets.For a three-dimensional description, the simplest basis set is obtained by a tensor prod-uct of one-dimensional basis functions. For a two resolution level description, the coarsedegrees of freedom are expanded by a single three dimensional function, φ i ,i ,i ( r ) , whilethe fine degrees of freedom can be expressed by adding another seven basis functions, φ νj ,j ,j ( r ) , which include tensor products with one-dimensional wavelet functions. Thus,Figure 4. Daubechies scaling function φ ( x ) and wavelet ϕ ( x ) of order 16.the Kohn-Sham wave function ψ ( r ) is of the form ψ ( r ) = X i ,i ,i c i ,i ,i φ i ,i ,i ( r ) + X j ,j ,j X ν c νj ,j ,j φ νj ,j ,j ( r ) . (52)The sum over i ,i ,i runs over all the grid points contained in the low-resolution regionsand the sum over j ,j ,j runs over all the points contained in the (generally smaller)high resolution regions. Each wave function is then described by a set of coefficients { c νj ,j ,j } , ν = 0 , ..., . Only the nonzero scaling function and wavelet coefficients arestored. The data is thus compressed. The basis set being orthogonal, several operationssuch as scalar products among different orbitals and between orbitals and the projectors ofthe nonlocal pseudopotential can be directly carried out in this compressed form. The matrix elements of the kinetic energy operator among the basis functions of our mixedrepresentation (i.e., scaling functions with scaling functions, scaling function with waveletsand wavelets with wavelets) can be calculated analytically [77]. For simplicity, let us il-lustrate the application of the kinetic energy operator onto a wavefunction ψ that is onlyexpressed in terms of scaling functions. ψ ( x, y, z ) = X i ,i ,i s i ,i ,i φ ( x/h − i ) φ ( y/h − i ) φ ( z/h − i ) avelets for DFT and Post-DFT Calculations 23The result of the application of the kinetic energy operator on this wavefunction, projectedto the original scaling function space, has the expansion coefficients ˆ s i ,i ,i = − h Z φ ( x/h − i ) φ ( y/h − i ) φ ( z/h − i ) ××∇ ψ ( x, y, z )dxdydz . Analytically the coefficients s i ,i ,i and ˆ s i ,i ,i are related by a convolution ˆ s i ,i ,i = 12 X j ,j ,j K i − j ,i − j ,i − j s j ,j ,j (53)where K i ,i ,i = T i T i T i , (54)where the coefficients T i can be calculated analytically via an eigenvalue equation: T i = Z φ ( x ) ∂ ∂x φ ( x − i ) dx = X ν,µ h ν h µ Z φ (2 x − ν ) ∂ ∂x φ (2 x − i − µ ) dx = X ν,µ h ν h µ − Z φ ( y − ν ) ∂ ∂y φ ( y − i − µ ) dy = X ν,µ h ν h µ Z φ ( y ) ∂ ∂y l φ ( y − i − µ + ν ) dy = X ν,µ h ν h µ T i − ν + µ Using the refinement equation Eq. (51), the values of the T i can be calculated analytically,from a suitable eigenvector of a matrix derived from the wavelet filters [77]. For this reasonthe expression of the kinetic energy operator is exact in a given Daubechies basis.Since the 3-dimensional kinetic energy filter K i ,i ,i is a product of three one-dimensional filters Eq. (54) the convolution in Eq. (53) can be evaluated with N N N L operations for a three-dimensional grid of N N N grid points. L is the length of theone-dimensional filter which is 29 for our Daubechies family. The kinetic energy can thusbe evaluated with linear scaling with respect to the number of nonvanishing expansion co-efficients of the wavefunction. This statement remains true for a mixed scaling function-wavelet basis where we have both nonvanishing s and d coefficients and for the case wherethe low and high resolution regions cover only parts of the cube of N N N grid points. In spite of the striking advantages of Daubechies wavelets the initial exploration of thisbasis set [78] did not lead to any algorithm that would be useful for practical electronicstructure calculations. This was due to the fact that an accurate evaluation of the localpotential energy is difficult in a Daubechies wavelet basis.4 Natarajan, Genovese, Casida, and DeutschBy definition, the local potential v ( r ) can be easily known on the nodes of the uniformgrid of the simulation box. Approximating a potential energy matrix element v i,j,k ; i ′ ,j ′ ,k ′ v i,j,k ; i ′ ,j ′ ,k ′ = Z d r φ i ′ , j ′ , k ′ ( r )v( r ) φ i , j , k ( r ) by v i,j,k ; i ′ ,j ′ ,k ′ ≈ X l,m,n φ i ′ ,j ′ ,k ′ ( r l,m,n ) v ( r l,m,n ) φ i,j,k ( r l,m,n ) gives an extremely slow convergence rate with respect to the number of grid points used toapproximate the integral because a single scaling function is not very smooth, i.e., it has arather low number of continuous derivatives. A. Neelov and S. Goedecker [79] have shownthat one should not try to approximate a single matrix element as accurately as possiblebut that one should try instead to approximate directly the expectation value of the localpotential. The reason for this strategy is that the wavefunction expressed in the Daubechybasis is smoother than a single Daubechies basis function. A single Daubechies scalingfunction of order 16 (i.e., the corresponding wavelet has 8 vanishing moments) has only 2continuous derivatives. More precisely its index of H ¨older continuity is about 2.7 and theSobolev space regularity with respect to p = 2 is about 2.91 [80]. A single Daubechiesscaling function of order 16 has only 4 continuous derivatives. By suitable linear com-binations of Daubechies 16 one can however exactly represent polynomials up to degree7, i.e., functions that have 7 non-vanishing continuous derivatives. The discontinuities getthus canceled by taking suitable linear combinations. Since we use pseudopotentials, ourexact wavefunctions are analytic and can locally be represented by a Taylor series. Weare thus approximating functions that are approximately polynomials of order 7 and thediscontinuities nearly cancel.Instead of calculating the exact matrix elements we therefore use matrix elements withrespect to a smoothed version ˜ φ of the Daubechies scaling functions. v i,j,k ; i ′ ,j ′ ,k ′ ≈ X l,m,n ˜ φ i ′ ,j ′ ,k ′ ( r l,m,n ) v ( r l,m,n ) ˜ φ i,j,k ( r l,m,n )= X l,m,n ˜ φ , , ( r l − i ′ ,m − j ′ ,n − k ′ ) V ( r l,m,n ) ˜ φ , , ( r l − i,m − j,n − k ) (55)where the smoothed wave function is defined by ˜ φ , , ( r l,m,n ) = ω l ω m ω n and ω l is the “magic filter”. Even though Eq. (55) is not a particulary good approximationfor a single matrix element it gives an excellent approximation for the expectation values ofthe local potential energy Z dx Z dy Z dzψ ( x, y, z ) v ( x, y, z ) ψ ( x, y, z ) and also for matrix elements between different wavefunctions Z dx Z dy Z dzψ i ( x, y, z ) v ( x, y, z ) ψ j ( x, y, z ) avelets for DFT and Post-DFT Calculations 25in case they are needed. Because of this remarkable achievement of the filter ω we call itthe magic filter.Following the same guidelines as the kinetic energy filters, the smoothed real spacevalues ˜ ψ i,j,k of a wavefunction ψ are calculated by performing a product of three one-dimensional convolutions with the magic filters along the x , y and z directions. For thescaling function part of the wavefunction the corresponding formula is ˜ ψ i ,i ,i = X j ,j ,j s j ,j ,j v (1) i − j v (1) i − j v (1) i − j where v (1) i is the filter that maps a scaling function on a double resolution grid. Similarconvolutions are needed for the wavelet part. The calculation is thus similar to the treatmentof the Laplacian in the kinetic energy.Once we have calculated ˜ ψ i,j,k the approximate expectation value ǫ V of the local po-tential v for a wavefunction ψ is obtained by simple summation on the double resolutionreal space grid: ǫ v = X j ,j ,j ˜ ψ j ,j ,j v j ,j ,j ˜ ψ j ,j ,j The energy contributions from the non-local pseudopotential have for each angular moment l the form X i,j h ψ | p i i h ij h p j | ψ i where | p i i is a pseudopotential projector. Once applying the Hamiltonian operator, theapplication of one projector on the wavefunctions requires the calculation of | ψ i → | ψ i + X i,j | p i i h ij h p j | ψ i . If we use for the projectors the representation of Eq. (52) (i.e., the same as for the wave-functions) both operations are trivial to perform. Because of the orthogonality of the basisset we just have to calculate scalar products among the coefficient vectors and to update thewavefunctions. The scaling function and wavelet expansion coefficients for the projectorsare given by [81] Z p ( r ) φ i ,i ,i ( r )d r , Z p( r ) ϕ ν i , i , i ( r )d r . (56)The GTH-HGH pseudopotentials [82, 83] have projectors which are written in termsof gaussians times polynomials. This form of projectors is particularly convenient to beexpanded in the Daubechies basis. In other terms, since the general form of the projector is h r | p i = e − cr x ℓ x y ℓ y z ℓ z , Z h r | p i φ i ,i ,i ( r )d r = W i ( c, ℓ x ) W i ( c, ℓ y ) W i ( c, ℓ x ) , (57) W j ( c, ℓ ) = Z + ∞−∞ e − ct t ℓ φ ( t/h − j )dt (58)The one-dimensional integrals are calculated in the following way. We first calculatethe scaling function expansion coefficients for scaling functions on a one-dimensional gridthat is 16 times denser. The integration on this dense grid is done by the well-known quadra-ture introduced in [84], that coincides with the magic filter [79]. This integration schemebased on the magic filter has a convergence rate of h and we gain therefore a factor of in accuracy by going to a denser grid. This means that the expansion coefficients arefor reasonable grid spacings h accurate to machine precision. After having obtained theexpansion coefficients with respect to the fine scaling functions we obtain the expansioncoefficients with respect to the scaling functions and wavelets on the required resolutionlevel by one-dimensional fast wavelet transformations. No accuracy is lost in the wavelettransforms and our representation of the projectors is therefore typically accurate to nearlymachine precision. In order to treat with the same advantages other pseudopotentials whichare not given under the form of gaussians it would be necessary to approximate them by asmall number of gaussians. Solving the Poisson equation for an arbitrary charge distribution is a non-trivial task, andis of major importance in many field of science, especially in the field of computationalchemistry. A huge effort has been put into making efficient Poisson solvers, and usualreal-space approaches includes finite difference (FD) and finite element (FE) methods. FDis a grid-based method, which is solving the equations iteratively on a discrete grid ofpointvalues, while FE is expanding the solution in a basis set, usually by dividing spaceinto cubic cells and allocate a polynomial basis to each cell.It is well-known fact that the electronic density in molecular systems is rapidly varyingin the vicinity of the atomic nuclei, and a usual problem with real-space methods is that anaccurate treatment of the system requires high resolution of gridpoints (FD) or cells (FE) inthe nuclear regions. Keeping this high resolution uniformly throughout space would yieldunnecessary high accuracy in the interatomic regions, and the solution of the Poisson equa-tion for molecular systems is demanding a multiresolution framework in order to achievenumerical efficiency. This chapter is concerned with a way of doing DFT and TD-DFT cal-culations, one where the multiresolution character is inherent in the theory, namely usingwavelet bases.In order to evaluate the Hartree potential, we need to rewrite the standard Poisson equa-tion to an integral form. The equation, in its differential form, is given as ∇ v ( x ) = 4 πρ ( x ) , (59)avelets for DFT and Post-DFT Calculations 27where ρ ( x ) is the known (charge) distribution, and v ( x ) is the unknown (electrostatic)potential. It is a standard textbook procedure to show that the solution can be written as theintegral v ( x ) = Z G ( x , y ) ρ ( y ) d y , (60)where G ( x , y ) is the Green’s function which is the solution to the fundamental equationwith homogeneous (Dirichlet) boundary conditions ∇ G ( x , y ) = δ ( x − y ) G ( x , y ) = 0 , x ∈ boundary (61)This equation can be solved analytically and the Green’s function is given (in three dimen-sions) simply as G ( x , y ) = 1 || x , y || , (62)This is the well-known potential arising from a point charge located in the position y , whichis exactly what Eq. (61) describes. The Green’s function kernel as it is given in Eq. (62) is not separable in the cartesian co-ordinates. However, since we are working with finite precision we can get by with an approximate kernel as long as the error introduced with this approximation is less than ouroverall accuracy criterion. If we are able to obtain such a numerical separation of the kernel,the operator can be applied in one direction at the time, allowing us to use the expressionsderived above for one-dimentional integral operators to solve the three-dimensional Pois-son equation. This is of great importance because it reduces the scaling behavior to becomelinear in the dimension of the system.The Poisson kernel can be made separable by expanding it as a sum of Gaussian func-tions, specifically r ≃ M ǫ X k =1 ω k e − p k r . (63)where ω k and p k are parameters that needs to be determined, and the number of terms M ǫ , called the separation rank, depends on the accuracy requirement and on what intervalthis expression needs to be valid. Details of how to obtain this expression can be found in[71, 72], and will not be treated here, but it should be mentioned that the separation rank isusually in the order of 100, e.g, it requires M ǫ = 89 to reproduce r on the interval [10 − , in three dimensions with error less than ǫ = 10 − .Finally, figure 5 summarize this complete section into a flow-chart type diagram. Thiskind of explanation is necessary for begineers because there are different functions used forthe different operations in B IG DFT. As one can see from the figure, The KS wavefunctions | ψ i are expressed in terms of Daubechies wavelets and the projection of Hamiltonian V nl | ψ i and of pseudopotential operators | Hψ i also expressed using Daubechies wavelets. The restof the operations such as kinetic energy |∇ ψ i , potential energy operator V ( x ) ψ ( x ) , and8 Natarajan, Genovese, Casida, and Deutschthe local densities ρ ( x ) are all expressed using interpolating scaling functions, in whichthe Hartree V H ( x ) , local part potential energy V loc ( x ) and xc operations V xc ( x ) were per-formed in real space. The interconnecting lines between different operations represents thetransformation between Daubechies wavelets-to-ISFs or the transformation of real space-to-fourier space representation.Figure 5. Operations performed in B IG DFT B IG DFT and TD-DFT
We want to solve Casida’s equation [48], (cid:20)(cid:18) A ( ω ) B ( ω ) B ∗ ( ω ) A ∗ ( ω ) (cid:19) − ω (cid:18) − (cid:19)(cid:21) (cid:18) XY (cid:19) = 0 , (64)where X and Y represents the pseudo eigenvectors; the matrices A and B are defined as A aiσ,bjτ = δ ab δ ij δ στ ( ǫ a − ǫ i ) + K aiσ,bjτ ( ω ) , (65)and, B aiσ,bjτ = K aiσ,jbτ ( ω ) , (66)in which the integral form of the coupling matrix K is given by, K pqσ,rsτ = Z Z Ψ ∗ pσ ( ~r )Ψ qσ ( ~r ) (cid:20) | ~r − ~r ′ | + ∂ E xc [ ρ ] ∂ρ σ ( ~r ) ∂ρ τ ( ~r ′ ) (cid:21) Ψ rτ ( ~r ′ )Ψ ∗ sτ ( ~r ′ ) d~rd~r ′ . (67)avelets for DFT and Post-DFT Calculations 29The universal adiabatic approximation is applied to Eq. (67) to remove the frequency de-pendence of the kernel.The electronic transitions occur with an infinitesimal perturbation obtains the abovedescribed non-Hermitian eigenvalue Eq. (64). Where the response is due to a real spinindependent external perturbation, and the actual response is described as the real densityresponse. However, an unitary transformation is necessary to convert Eq. (64) into the realeiganvalue problem. In Eq. (64), all occupied-occupied and virtual-virtual element contri-butions are zero whereas only the elements that are from virtual-occupied and occupied-virtual parts are taken into account. Moreover if we only restricted to virtual-occupiedelements and neglecting the occupied-virtual elements of Eq. (64) leads to a Hermitianeigenvalue equation of the dimension one-half of that TD-DFT working equation is said tobe Tamm-Dancoff approximation (TDA) and it is written as, AX = ω X , (68)where A is as same as in Eq. (42). The matrix A is just restricted to number of singleexcitations.However, the explicit form of Eq. (68) is, Ω( ω ) ~F I = ω ~F I , (69)where Ω iaσ,jbτ = δ ia δ jb δ στ ( ǫ aσ − ǫ iσ ) + (70) p ( ǫ aσ − ǫ iσ ) K iaσ,jbτ p ( ǫ aσ − ǫ iσ ) , where ǫ iσ − ǫ aσ is the energy eigenvalue differences of i th and a th states. Solving Eqs. (69)yields TD-DFT excitation energies ω and ~F I ’s are the corresponding oscillator strengthswhich are defined from the transition dipole moments. We are now in a position to understand the construction of the coupling matrix Eq. (67) inour implementation of TD-DFT in B IG DFT, which we split into the Hartree andxc parts, K aiσ,bjτ = K Haiσ,bjτ + K xcajσ,bjτ . (71)Instead of calculating the Hartree part of coupling matrix directly as, K Haiσ,bjτ = Z Z ψ ∗ aσ ( r ) ψ iσ ( r ) 1 | r − r ′ | ψ bτ ( r ′ ) ψ ∗ jτ ( r ′ ) d r d r ′ , (72)we express the coupling matrix element as, K Haiσ,bjτ = Z ψ ∗ aσ ( r ) ψ iσ ( r ) v bjτ ( r ) d r , (73)where, v aiσ ( r ) = Z ρ aiσ ( r ) | r − r ′ | d r ′ , (74)0 Natarajan, Genovese, Casida, and Deutschand, ρ aiσ ( r ) = ψ ∗ aσ ( r ) ψ iσ ( r ) . (75)The advantage of doing this is that, although ρ aiσ and v aiσ are neither real physical chargedensities nor real physical potentials, they still satisfy the Poisson equation, ∇ v aiσ ( r ) = − πρ aiσ ( r ) , (76)and we can make use of whichever of the efficient wavelet-based Poisson solvers alreadyavailable in B IG DFT, is appropriate for the boundary conditions of our physical problem.Once the solution of Poisson’s equation, v aiσ ( r ) , is known, we can then calculate theHartree part of the kernel according to Eq. (73). Inclusion of the xc kernel is accomplishedby evaluating, K aiσ,bjτ = Z M aiσ ( r ) ρ bjτ ( r ) d r , (77)where, M aiσ ( r ) = v aiσ ( r ) + Z ρ aiσ ( r ′ ) f σ,τxc ( r , r ′ ) d r ′ . (78)We note that f σ,τxc ( r , r ′ ) = f σ,τxc ( r , r ′ ) δ ( r − r ′ ) for the LDA, so that no integral need actuallybe carried out in evaluating M aiσ ( r ) . The integral in Eq. (77) is, of course, carried outnumerically in practice as a discrete summation.
8. Results
We now wish to illustrate a bit how wavelet calculations work in the B IG DFT pro-gram. Comparison will be made against results obtained with the GTO-based program DE M ON K . This work is very similar to our previous work reporting the first implemen-tation of wavelet-based TD-DFT with illustration for N and application to the absorptionspectrum of a medium-sized organic molecule of potential biomedical use as a fluores-cent probe [45]. Here however we will present new B IG DFT results for a different smallmolecule, namely carbon monoxide. Though CO is roughly isoelectronic with N , CO hasthe interesting feature of having a low-lying bright state in its absorption spectrum. Calculations were carried out with DE M ON K and B IG DFT with the LDA-optimized bondlength of 1.129 ˚A. DE M ON KDE M ON K resembles a typical GTO-based quantum chemistry program in that all the in-tegrals other than the xc-integrals, can be evaluated analytically. In particular, DE M ON K has the important advantage that it accepts the popular GTO basis sets common in quan-tum chemistry and so can benefit from the experience in basis set construction of a largecommunity built up over the past 50 years or so. In the following, we have chosen to usethe well-known correlation-consistent basis sets for this study [85, 86]. (Note, however,avelets for DFT and Post-DFT Calculations 31that the correlation-consistent basis sets used in DE M ON K lack f and g functions but areotherwise exactly the same as the usual ones.) The advantage of using these particular basissets is that there is a clear hierarchy as to quality.An exception to the rule that integrals are evaluated analytically in DE M ON K are thexc-integrals (for the xc-energy, xc-potential, and xc-kernel) which are evaluated numeri-cally over a Becke atom-centered grid. This is important because the relative simplicity ofevaluating integrals over a grid has allowed the rapid implemenation of new functionals asthey were introduced. We made use of the fine fixed grid in our calculations.As described so far, DE M ON K should have O ( N ) scaling because of the need toevaluate 4-center integrals. Instead DE M ON K uses a second atom-centered auxiliary GTObasis to expand the charge density. This allows the the elimination of all 4-center integralsso that only 3-center integrals remain for a formal O ( N ) scaling. In practice, integralprescreening leads to O ( N M ) scaling where M is typically between 2 and 3. We made useof the A3 auxiliary basis set from the DE M ON K automated auxiliary basis set library.All calculations were performed using standard DE M ON K default criteria. The im-plementation of TD-DFT in DE M ON K is described in Ref. [87]. (The charge densityconservation constraint is no longer used in DE M ON K TD-DFT calculations.) Althoughfull TD-LDA calculations are possible with DE M ON K , the TD-LDA calculations reportedhere all made use of the TDA. B IG DFT (a) H O in a box (b) H O in a box showing fine gridresolution (c) H O in a box showing coarsegrid resolution
Figure 6. Adaptive grid in B IG DFT (a), (b) and (c)The main thing to vary in B IG DFT is the grid which is of more profound importancethan in DE M ON K because it is the grid which supports the wavelets. Figures 6(a), 6(b),and 6(c) give an idea of what the grid looks like for the small familiar molecule of water.Conceptually the molecule is in a very large box (Fig. 6(a).) A fine grid is placed in theregions of high electron density around the molecule (Fig. 6(b).) A coarse grid is used ina larger region where the electron density varies more slowly (Fig. 6(c).) The B IG DFTgrid is characterized by the triple h g /crmult/frmult. The first number in the triple ( h g ) isa real number which specifies the nodes of the grid in atomic units. The second number(the integer-valued crmult) is the coarse grid multiplier. And the third number (the integer-valued frmult) is the fine grid multiplier. Two points must be clearly understood whenlooking at Figs. 6(a), 6(b), and 6(c). The first is that, while the box may determine the limits2 Natarajan, Genovese, Casida, and DeutschTable 1. Basis set dependence of the HOMO and LUMO energies and of the HOMO-LUMO gap (eV) calculated using DE M ON K . Basis Set − ǫ HOMO − ǫ LUMO ∆ ǫ HOMO − LUMO
STO-3G -5.5350 1.2428 4.2922DZVP -8.9271 -2.0942 6.8329TZVP -9.0287 -2.1902 6.8385CC-PVDZ -8.6729 -1.7823 6.8906CC-PVTZ -9.0419 -2.1195 6.9224CC-PVQZ -9.0944 -2.1971 6.8973CC-PV5Z -9.1169 -2.2400 6.8769CC-PCVDZ -8.6905 -1.7922 6.8983CC-PCVQZ -9.0957 -2.1988 6.8969CC-PCVTZ -9.0371 -2.1165 6.9206CC-PCV5Z -9.1172 -2.2401 6.8771AUG-CC-PVDZ -9.0910 -2.2345 6.8565AUG-CC-PVQZ -9.1286 -2.2567 6.8719AUG-CC-PVTZ -9.1306 -2.2535 6.8771AUG-CC-PV5Z -9.1289 -2.2606 6.8683AUG-CC-PCVDZ -9.0987 -2.2371 6.8616AUG-CC-PCVTZ -9.1316 -2.2554 6.5776AUG-CC-PCVQZ -9.1293 -2.2574 6.8719AUG-CC-PCV5Z -9.1291 -2.2607 6.8684 of the grid, the grid does not have the shape of the box and there are no basis functions wherethere are no grid points. This means that we are not dealing with box boundary conditions,but rather with effective boundary conditions which reflect the shape of the molecule. Theother point which is not brought out by our explanation is that the B IG DFT grid is adaptivein the sense that additional fine grid points are added during the calculation as they areneeded to maintain and improve numerical precision.The implementation of TD-DFT in B IG DFT is described in Ref. [45].
Possibly the most remarkable property of wavelets is how rapidly they converge to thebasis set limit. Let us illustrate this by comparing highest-occupied molecular orbital(HOMO) and lowest-unoccupied molecular orbital (LUMO) energies calculated with DE -M ON K and B IG DFT. The difference of these two energies is the HOMO-LUMO gap, ∆ ǫ HOMO − LUMO .Consider first how DE M ON K calculations of ∆ ǫ HOMO − LUMO , evolve as the basisset is improved (Table 1.) Convergence to the true HOMO-LUMO LDA gap is expectedavelets for DFT and Post-DFT Calculations 33Table 2. Basis set dependence of the HOMO and LUMO energies and of the HOMO-LUMO gap (eV) calculated using B IG DFT. h ga /m b /n c − ǫ HOMO − ǫ LUMO ∆ ǫ HOMO − LUMO a Grid spacing of the cartesian grid in atomic units. b Coarse grid multiplier (crmult). c Fine grid multiplier (frmult). with systematic improvement within the series: (i) double zeta plus valence polariza-tion (DZVP) → triple zeta plus valence polarization (TZVP), (ii) augmented correlation-consistent double zeta plus polarization plus diffuse on all atoms (AUG-CC-PCVDZ) → AUG-CC-PCVTZ (triple zeta) → AUG-CC-PCVQZ (quadruple zeta) → AUG-CC- PCV5Z(quintuple zeta), (iii) augmented correlation-consistent valence double zeta plus polariza-tion plus diffuse (AUG-CC-PVDZ) → AUG-CC-PVTZ → AUG-CC-PVQZ → AUG-CC-PV5Z, (iv) correlation-consistent double zeta plus polarization plus tight core (CC-PCVDZ) → CC-PCVTZ → CC-PCVQZ → CC-PCV5Z, and (v) correlation-consistent valence dou-ble zeta plus polarization on all atoms (CC-PVDZ) → CC-PVTZ → CC-PVQZ → CC-PV5Z. There is a clear tendency in the correlation-consistent basis sets to tend towardsvalues of -9.13 eV for the HOMO energy, -2.26 eV for the LUMO energy, and 6.87 eV for ∆ ǫ HOMO − LUMO , with adequate convergance achieved with the AUG-CC-PVQZ basis set.Now let us turn to BigDFT (Table 2). Calculations were done for several different grids,including the high-resolution combination 0.3/8/8 and the low-resolution combination of0.4/6/8. Remarkably, except for the very lowest quality grid 0.4/6/8, there is essentiallyno difference between results obtained with the two grids (and even the 0.4/6/8 grid givesnearly converged results.) The results are also quite close to, but not identical to thoseobtained with the DE M ON K program. The reason for the small differences between theconverged results obtained with the two programs is more difficult to trace as it might bedue to the auxiliary basis approximation in DE M ON K or to the use of pseudopotentials inB IG DFT or perhaps to still other program differences. The important point is that differ-ences are remarkably small.
Orbital energy differences provide a first estimate for excitation energies. In this case, wewould expect to see the HOMO → LUMO excitation at ∆ ǫ HOMO − LUMO ≈ . eV (6.87eV for DE M ON K and 6.90 eV for B IG DFT.) A better estimate is provided by the two-orbital two-electron model (TOTEM) [48, 88, 89, 90] for the singlet (S) and triplet (T)4 Natarajan, Genovese, Casida, and DeutschTable 3. Comparison of lowest excitation energies of CO (in eV) calculated using B IG DFTand DE M ON K and with experiment. State B IG DFT a DE M ON K b Experiment c Σ − ∆ Π Σ + Π a Present work (TD-LDA/TDA) using AUG-CC-PCQZ basis set. b Present work (TD-LDA/TDA) using 0.3/8/8 grid. c Taken from Ref. [91]. transition from orbital i to orbital a , ~ ω Ti → a = ∆ ǫ i → a + ( ia | f α,αxc − f α,βxc | ai ) ~ ω Si → a = ∆ ǫ i → a + ( ia | f H + f α,αxc + f α,βxc | ai ) , (79)where ∆ ǫ i → a = ǫ a − ǫ i . (80)The TOTEM model often works surprisingly well for small molecules because, unlikethe Hartree-Fock approximation which is better adapted to describe electron ionization andattachment, pure DFT Kohn-Sham orbitals are preprepared to describe excitation energiesin the sense that the occupied and unoccupied orbitals see the same potential, thus minimiz-ing orbital relaxation effects. Inspection of the sizes and signs of the integrals in Eq. (79)indicates that we should expect, ~ ω Ti → a ≤ ∆ ǫ i → a ≤ ~ ω Si → a . (81)These is confirmed in Table 3 where the Π and Π excitations are, respectively, thetriplet and singlet states corresponding to the HOMO → LUMO transition.Assuming that f α,αxc dominates over f α,βxc , we may even go a bit further to estimate ( ia | f H | ai ) and ( ia | f α,αxc | ai ) (Fig. 7.) The calculations are show in Table 4. Comparison of ( ia | f α,αxc | ai ) (1) and ( ia | f α,αxc | ai ) (2) provides an indication of the quality of the approxima-tion of neglecting the ( ia | f α,βxc | ai ) integral which in this case appears to be excellent. The ( ia | f H | ai ) integrals calculated with the two programs are reasonably close. Interestinglythe ( ia | f α,αxc | ai ) disagree by about 0.4 eV which, though small, is not negligible.Let us now examine the issue of the collapse of the continuum. In Ref. [92], it wasshown that the TD-DFT ionization continuum begins at − ǫ HOMO . In exact Kohn-ShamDFT, this should be the ionization potential. However typical approximate density func-tionals underbind electrons and so lead to an artificially-early on-set of the TD-DFT ion-ization continuum. This is first illustrated using the DE M ON K program and different basissets. Indeed Fig. 8 shows that the states above − ǫ HOMO tend to collapse towards − ǫ HOMO rather than converging as they should. This is simply because we are trying to describe acontinuum which should not be there with a finite basis set. Also seen in the figure is aavelets for DFT and Post-DFT Calculations 35Figure 7. Estimation of integrals within the TOTEM model.Table 4. Estimations of integrals (in eV) within the TOTEM.
Program DE M ON K B IG DFTInput Data Π ∆ ǫ i → a Π ( ia | f H | ai ) ( ia | f α,αxc | ai ) (1) -0.82 -0.43 ( ia | f α,αxc | ai ) (2) -0.83 -0.44 slight splitting of the Π excitation energy. This small effect is due to the fact that the gridused to calculate xc-integrals in DE M ON K has only roughly the symmetry of the molecule.Now let us turn to B IG DFT calculations. Figure 9 shows a similar collapse of the con-tinuum as the fineness of the grid increases. Interestingly there is no evidence of symmetrybreaking of the doubly-degenerate Π state. Carbon monoxide is very unusual for small molecules in that absolute oscillator strengthshave been well studied [93] over a significant energy range and the A Π ( Π in Table 3)is bright and has an accurately determined oscillator strength. See Fig. 2 of Ref. [93] (aswell as other references in the same paper) for a graph of measured absolute optical oscilla-tor strengths against absorption energy in eV. Table 5 reports our calculated TD-LDA/TDA6 Natarajan, Genovese, Casida, and Deutsch STO-3G TZVP CC-PVTZ CC-PCVTZ AUG-CC-PVTZ AUG-CC-PCVTZ E x c it a ti on e n e r g i e s ( e V ) - ε HOMO
Figure 8. Singlet and triplet excitation energies for CO calculated using DE M ON K E x c it a ti on e n e r g i e s ( e V ) - ε HOMO
Figure 9. Singlet and triplet excitation energies for CO calculated using B IG DFToscillator strengths. As the TDA violates the Thomas-Reiche-Kuhn (TRK) f -sum rule[48] it should only be used very cautiously to estimate oscillator strengths. Neverthelessthe DE M ON K value of f = 0 . is in good agreement with the experimental value of f = 0 . . As shown in Ref. [91], full TDLDA calculations with asymptotically cor-rected potentials give smaller oscillator strengths (0.136 for TD-LDA/LB94 and 0.156 forTD-LDA/AC-LDA calculations[91].) (Co¨ıncidently our own DE M ON K full TD-LDA cal-cuations without asymptotic corrections give a degeneracyc-weighted oscillator strength of0.1752 (bang on the experimental value) but an excitation energy of 8.19 eV.) Since os-avelets for DFT and Post-DFT Calculations 37Table 5. Comparison of experimental A Π energies (eV) and oscillator strengths with TD-LDA/TDA experimental A Π energies (eV) and degeneracy-weighted oscillator strengths(unitless.) DE M ON K B IG DFT Experiment a ~ ω S f a See Table VIII of Ref. [91].cillator strengths are quite sensitive to configuration mixing with nearby states, the factthat the B IG DFT oscillator strength is larger than the DE M ON K oscillator strength maybe due to the relatively small energy separation between the B IG DFT A Π state and theartificially-low TD-LDA ionization continuum.
9. Conclusion
In this chapter we have tried to give an informative elementary review of a subject largelyunfamiliar to most theoretical chemists and physicists. Wavelets, once an obscure ripple atthe exterior of engineering applications, grew to become a regular tsumani in engineeringcircles in the 1990s as the similarity to and superiority over Fourier transform methods formultiresolution problems with arbitrary boundary conditions became increasingly recog-nized. Though the first applications of wavelet theory to solving the Schr¨odinger equationmay be traced back to the mid-1990s [94, 95, 96], the theory is still not well known amongquantum mechanicians. Here we have tried to remedy this aberrant situation by trying to“make some waves about wavelets for wave functions.”In particular we have reviewed the theory behind the wavelet code B IG DFT for ground-state DFT and our recent implementation of wavelet-based TD-DFT in B IG DFT. Rapidprogress is being made towards making B IG DFT a high performance computing order- N code for applications to large systems. Right now applications to 400 or 500 atomsare routine for ground-state calculations with B IG DFT. Our implementation of TD-DFTin B IG DFT is by comparison only a rudementary beginning, but it shows that the basicmethod is viable and we are confident that there are no insurmountable obstacles to makinghigh performance computing order- N wavelet-based TD-DFT code for large systems. Acknowledgments
B. N. would like to acknowledge a scholarship from the
Foundation Nanoscience . Thiswork has been carried out in the context of the French Rhˆone-Alpes
R´eseau th´ematiquede recherche avanc´ee (RTRA): Nanosciences aux limites de la nano´electronique and theRhˆone-Alpes Associated Node of the European Theoretical Spectroscopy Facility (ETSF).8 Natarajan, Genovese, Casida, and Deutsch
AppendicesA List of Abbreviations
For the readers convenience we give a list of the abbreviations used in this chapter in alpha-betical order: AA Adiabatic approximation.
APW
Augmented plane wave. CC Coupled cluster. CI Configuration interaction.
DFT
Density-functional theory. DM Density matrix. FD Finite difference. FE Finite element.
FEM
Finite element method.
GTH-HGH
Goedecker-Teter-Hutter/Hartwigsen-Goedecker-Hutter.
GTO
Gaussian-type orbitals. H Hartree. HF Hartree-Fock.
HOMO
Highest-occupied molecular orbital.
ISF
Interpolating scaling function.
LCAO
Linear combination of atomic orbitals.
LDA
Local density approximation.
LMTO
Linear muffin tin orbital. LR Linear-response.
LR-TD-DFT
Linear-response time-dependent density-functional theory.
LR-TD-LDA
Linear-response time-dependent local density approximation.
LUMO
Lowest-unoccupied molecular orbital. NS Non-standard.
MBPT
Many-body perturbation theory.avelets for DFT and Post-DFT Calculations 39
MRA
Multiresolution analysis. S Singlet.
SOS
Sum-over-states.
STO
Slater-type orbital T Triplet. TD Time-dependent.
TDA
Tamm-Dancoff approximation.
TD-DFT
Time-dependent density-functional theory.
TD-LDA
Time-dependent local density approximation.
TRK
Thomas-Reiche-Kuhn. xc Exchange-correlation.
References [1] T. H. Dunning,
Gaussian basis sets for use in correlated molecular calculations.I. The atoms boron through neon and hydrogen , J. Chem. Phys. , 1007 (1989).[2] R. A. Kendall, T. H. Dunning, and R. J. Harrison, Electron affinities of the first-rowatoms revisited. Systematic basis sets and wave functions , J. Chem. Phys. ,6796 (1992).[3] D. E. Woon and T. H. Dunning, Gaussian basis sets for use in correlated molec-ular calculations. III. The atoms aluminum through argon , J. Chem. Phys. ,1358 (1993).[4] X. P. Li, R. W. Nunes, and D. Vanderbilt, Density-matrix electronic-structuremethod with linear system-size scaling , Phys. Rev. B , 10891 (1993).[5] M. S. Daw, Model for energetics of solids based on the density matrix , Phys.Rev. B , 10895 (1993).[6] S. Goedecker and L. Colombo, Efficient linear scaling algorithm for tight-bindingmolecular dynamics , Phys. Rev. Lett. , 122 (1994).[7] W. Kohn, Density functional and density matrix method scaling linearly withthe number of atoms , Phys. Rev. Lett. , 3168 (1996).[8] J. M. Millam and G. E. Scuseria, Linear scaling conjugate gradient density matrixsearch as an alternative to diagonalization for first principles electronic struc-ture calculations , J. Chem. Phys. , 5569 (1997).0 Natarajan, Genovese, Casida, and Deutsch[9] C. Ochsenfeld and M. Head-Gordon,
A reformulation of the coupled perturbedself-consistent field equations entirely within a local atomic orbital densitymatrix-based scheme , Chem. Phys. Lett. , 399 (1997).[10] B. Huron, J. P. Malrieu, and P. Rancurel,
Iterative perturbation calculationsof ground and excited state energies from multiconfigurational zeroth-orderwavefunctions , J. Chem. Phys. , 5745 (1973).[11] S. Evangelisti, J. P. Daudey, and J. P. Malrieu, Convergence of an improved CIPSIalgorithm , Chem. Phys. , 91 (1983).[12] M. Sch¨utz, G. Hetzer, and H. J. Werner, Low-order scaling local electron correla-tion methods. I. Linear scaling local MP2 , J. Chem. Phys. , 5691 (1999).[13] G. E. Scuseria and P. Y. Ayala,
Linear scaling coupled cluster and perturbationtheories in the atomic orbital basis , J. Chem. Phys. , 8330 (1999).[14] D. R. Hartree,
The Calculation of Atomic Structures , Wiley, New York, 1957.[15] V. Fock,
N ¨aherungsmethode zur L ¨osung des quantenmechanischenMehrk ¨orperproblems , Z. f¨ur physik , 126 (1930).[16] J. C. Slater, Note on Hartree’s method , Phys. Rev. , 210 (1930).[17] C. C. J. Roothan, New developments in molecular orbital theory , Rev. Mod. Phys. , 69 (1951).[18] L. H. Thomas, The calculation of atomic fields , Proc. Cambridge Phil. Soc. , 542(1926).[19] E. Fermi, Sulla quantizzazione del gas perfetto monoatomico , Rend. Lincei ,145 (1926).[20] E. Fermi, Zur Quantelung des idealen einatomigen Gases , Z. f¨ur physik , 902(1926).[21] J. Frenkel, Zur wellenmechanischen Theorie der metallischen Leitf ¨ahigkeit , Z.f¨ur physik , 819 (1928).[22] A. Sommerfeld, Zur Elektronentheorie der Metalle auf Grund der FermischenStatistik I. Teil: Allgemeines, Str ¨omungs- und Austrittsvorg ¨ange , Z. f¨ur physik , 1 (1928).[23] P. A. M. Dirac, On the theory of quantum mechanics , Proc. Roy. Soc. London A , 661 (1926).[24] E. Fermi,
Un metodo statistico per la determinazione di alcune priopriet `adell’atomo , Rend. Accad. Naz. Lincei , 602 (1927).[25] P. A. M. Dirac, Note on exchange phenomena in the thomas atom , Math. Proc.Cambr. Phil. Soc. , 376 (1930).avelets for DFT and Post-DFT Calculations 41[26] P. Hohenberg and W. Kohn, Inhomogeneous electron gas , Phys. Rev. , B864(1964).[27] W. Kohn and L. J. Sham,
Self-consistent equations including exchange-correlation effects , Phys. Rev. A , 1133 (1965).[28] R. M. Dreizler and E. K. U. Gross,
Density Functional Theory , Springer-Verlag,Berlin Heidelberg New York, 1990.[29] F. Herman, J. P. V. Dyke, and I. B. Ortenburger,
Improved statistical exchangeapproximation for inhomogeneous many-electron systems , Phys. Rev. Lett. ,807 (1969).[30] O. K. Andersen, Linear methods in band theory , Phys. Rev. B , 3060 (1975).[31] E. Wimmer, H. Krakauer, M. Weinert, and A. J. Freeman, Full-potential self-consistent linearized-augmented-plane-wave method for calculating the elec-tronic structure of molecules and surfaces: O molecule , Phys. Rev. B , 864(1981).[32] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Johannopoulos, Itera-tive minimization techniques for ab initio total-energy calculations: Moleculardynamics and conjugate gradients , Rev. Mod. Phys. , 1045 (1992).[33] J. Pipek and I. Varga, Statistical electron densities , Int. J. Quant. Chem. , 85(1997).[34] R. Poirier, R. Kari, and I. G. Csizmadia, Handbook of Gaussian Basis sets , Elsevier,Amsterdam, 1970.[35] J. A. Pople and D. L. Beveridge,
Approximate Molecular Orbital Theory , McGraw-Hill, New York, 1970.[36] S. Wei and M. Y. Chou,
Wavelets in self-consistent electronic structure calcula-tions , Phys. Rev. Lett. , 2650 (1996).[37] T. Kato, On the eigenfunctions of many particle systems in quantum mechan-ics , Comm. Pure Appl. Math , 151 (1957).[38] T. A. Arias, R. A. Lippert, and A. Edelman, Multiscale computation with interpo-lating wavelets , J. Comput. Phys. , 278 (1997).[39] J. C. Slater,
Quantum Theory of Atomic Structure , McGraw-Hill, New York, 1960.[40] C. F. Fischer,
Numerical solution of the Hartree-Fock equations , Can. J. Phys. , 1895 (1963).[41] C. F. Fischer, Average-energy-of-configuration Hartree-Fock results for theatoms helium to radon charlotte froese fischer , At. Data Nucl. Data Tables ,301 (1972).2 Natarajan, Genovese, Casida, and Deutsch[42] C. F. Fischer, Average-energy-of-configuration Hartree-Fock results for theatoms helium to radon , At. Data Nucl. Data Tables , 87 (1973).[43] J. B. Mann, SCF Hartree-Fock results for elements with two open shells andthe elements francium to nobelium , At. Data Nucl. Data Tables , 1 (1973).[44] D. Heinemann, B. Ficke, and D. Kolb, Accurate Hartree-Fock-Slater calculationson small diatomic molecules with the finite-element method , Chem. Phys. Lett. , 125 (1988).[45] B. Natarajan, L. Genovese, M. E. Casida, T. Deutsch, O. N. Burchak, C. Philouze,and M. Y. Balakirev,
Wavelet-based linear-response time-dependent density-functional theory , http://arxiv.org/abs/1108.3475.[46] E. Runge and E. K. U. Gross,
Density-functional theory for time-dependent sys-tems , Phys. Rev. Lett. , 997 (1984).[47] M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Excitation energies from time-dependent density-functional theory , Phys. Rev. Lett. , 1212 (1996).[48] M. E. Casida, Time-dependent density-functional response theory formolecules , in
Recent Advances in Density Functional Methods, Part I , edited byD. P. Chong, page 155, World Scientific, Singapore, 1995.[49] C. Jamorski, M. E. Casida, and D. R. Salahub,
Dynamic polarizabilities and ex-citation spectra from a molecular implementation of time-dependent density-functional response theory: N as a case study , J. Chem. Phys. , 5134 (1996).[50] R. Bauernschmitt and R. Ahlrichs, Treatment of electronic excitations within theadiabatic approximation of time-dependent density functional theory , Chem.Phys. Lett., , 454 (1996).[51] C. Bienia, S. Kumar, J. P. Singh, and K. Li,
The PARSEC benchmark suite: Char-acterization and architectural implications , Technical report, in Princeton Univer-sity, 2008.[52] A. Castro, M. A. L. Marques, H. Appel, M. Oliveira, C. A. Rozzi, X. Andrade,F. Lorenzen, E. K. U. Gross, and A. Rubio,
Octopus: A tool for the applicationof time-dependent density functional theory , Physica Status Solidi , 2465(2006).[53] L. Lehtovaara, V. Havu, and M. Puska,
All-electron density functitonal theory andtime-dependent density-functional theory with high-order finite elements , J.Chem. Phys. , 054103 (2009).[54] O. M. Nielsen,
Wavelets in scientific computing , PhD thesis, Department of Mathe-matical Modelling, Technical University of Denmark, 1998.[55] W. Dahmen, A. Cohen, and R. DeVore,
Adaptive wavelet schemes for ellipticoperator equations - convergence rates , Math. Comput. , 27 (2001).avelets for DFT and Post-DFT Calculations 43[56] T. K. Jensen, On adaptive wavelet-based methods for the Maxwell equations , PhDthesis, Department of Mathematics, Technical University of Denmark, 2003.[57] http://inac.cea.fr/L_Sim/BigDFT/ .[58] F. Keinert,
Wavelets and Multiwavelets , Chapman and Hall/CRC, 2003.[59] A. Haar,
Zur Theorie des orthogonalen Funktionensysteme , MathematischeAnnalen , 331 (1910).[60] B. K. Alpert, A class of bases in L for the sparse representation of integraloperators , SIAM J. Math. Anal. , 246 (1993).[61] I. Daubechies, Ten Lectures on Wavelets , SIAM, Philadelphia, 1992.[62] R. A. Lippert, T. A. Arias, and A. Edelman,
Multiscale computation with interpo-lating wavelets , J. Comput. Phys. , 278 (1998).[63] P. Pulay,
Convergence acceleration of iterative sequences. The case of SCFiteration , Chem. Phys. Lett. , 393 (1980).[64] E. R. Davidson, The iterative calculation of a few of the lowest eigenvalues andcorresponding eigenvectors of large real-symmetric matrices , J. Comput. Phys. , 87 (1975).[65] C. W. Murray, S. C. Racine, and E. R. Davidson, Improved algorithms for thelowest few eigenvalues and eigenvectors of large matrices , J. Comput. Phys. , 382 (1991).[66] B. Walker, A. M. Saitta, R. Gebauer, and S. Baroni,
Efficient approach to time-dependent density-functional perturbation theory for optical spectroscopy ,Phys. Rev. Lett. , 113001 (2006).[67] D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, Turbo charging time-dependentdensity-functional theory with Lanczos chains , J. Chem. Phys. , 154105(2008).[68] S. Baroni, R. Gebauer, O. B. Malcioˇglu, Y. Saad, P. Umari, and J. Xian,
Harness-ing molecular excited states with Lanczos chains , J. Phys. Condens. Matter ,074204 (2010).[69] O. B. Malcioˇglu, R. Gebauer, D. Rocca, and S. Baroni, TURBO
TDDFT – A codefor the simulation of molecular spectra using the Liouville-Lanczos approachto time-dependent density-functional perturbation theory , Comp. Phys. Comm. , 1744 (2011).[70] DE M ON K @G RENOBLE , the Grenoble development version of DE M ON K , AndreasM. K ¨oster, Patrizia Calaminici, Mark E. Casida, Roberto Flores-Morino, GeraldGeudtner, Annick Goursot, Thomas Heine, Andrei Ipatov, Florian Janetzko, Sergei4 Natarajan, Genovese, Casida, and DeutschPatchkovskii, J. Ulisis Reveles, Dennis R. Salahub, and Alberto Vela, The Interna-tional deMon Developers Community (Cinvestav-IPN, Mexico, 2006) plus some ad-ditional features by Mark E. Casida, Lo¨ıc Joubert Doriol, Andrei Ipatov, Miquel Huix-Rotllant, and Bhaarathi Natarajan (Grenoble, France, 2011).[71] L. Genovese, T. Deutsch, A. Neelov, S. Goedecker, and G. Beylkin,
Efficient solutionof poisson’s equation with free boundary conditions , J. Chem. Phys., , 074105(2006).[72] L. Genovese, T. Deutsch, and S. Goedecker,
Efficient and accurate three-dimensional poisson solver for surface problems , J. Chem. Phys. , 054704(2007).[73] L. Genovese et al.,
Daubechies wavelets as a basis set for density functionalpseudopotential calculations , J. Chem. Phys. , 014109 (2008).[74] L. Genovese, M. Ospici, T. Deutsch, J.-F. M´ehaut, A. Neelov, and S. Goedecker,
Den-sity functional theory calculation on many-cores hybrid cpu-gpu architectures ,J. Chem. Phys. , 034103 (2009).[75] T. Deutsch and L. Genovese,
Wavelets for electronic structure calculations ,Journ´ees des Neutrons , 33 (2011).[76] L. Genovese, B. Videau, M. Ospici, T. Deutsch, S. Goedecker, and J.-F. M´ehaut, Daubechies wavelets for high performance electronic structure calculations:The bigdft project , Comptes Rendus M´ecanique , 149 (2011).[77] G. Beylkin,
On the representation of operators in bases of compactly sup-ported wavelets , SIAM J. Numer. Anal. , 1716 (1992).[78] C. J. Tymczak and X.-Q. Wang, Orthonormal wavelet bases for quantum molec-ular dynamics , Phys. Rev. Lett. , 3654 (1997).[79] A. I. Neelov and S. Goedecker, An efficient numerical quadrature for the calcu-lation of the potential energy of wavefunctions expressed in the Daubechieswavelet basis , J. Comp. Phys. , 055501 (2006).[80] K.-A. Lau and Q. Sun,
Asymptotic regularity of Daubechies scaling functions ,Proc. Am. Math. Soc. , 1087 (2000).[81] S. Goedecker,
Wavelets and Their Application for the Solution of Partial DifferentialEquations , Presses Polytechniques Universitaires et Romandes, Lausanne, Switzer-land, 1998.[82] S. Goedecker, M. Teter, and J. Hutter,
Separable dual-space Gaussian pseudopo-tentials , Phys. Rev. B , 1703 (1996).[83] S. Goedecker, M. Teter, and J. Hutter, Relativistic separable dual-space Gaussianpseudopotentials from H to Rn , Phys. Rev. B , 3641 (1998).avelets for DFT and Post-DFT Calculations 45[84] B. R. Johnson, J. P. Modisette, P. J. Nordlander, and J. L. Kinsey, Quadrature inte-gration for orthogonal wavelet systems , J. Chem. Phys. , 8309 (1999).[85] D. Feller,
The role of databases in support of computational chemistry calcua-tions , J. Comp. Phys. , 1571 (1996).[86] K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li,and T. L. Windus, Basis set exchange: A community database for computa-tional sciences , J. Chem. Inf. Model. , 1045 (2007).[87] A. Ipatov, A. Fouqueau, C. P. del Valle, F. Cordova, M. E. Casida, A. M. K ¨oster,A. Vela, and C. J. Jamorski, Excitation energies from an auxiliary-function for-mulation of time-dependent density-functional response theory with chargeconservation constraint , J. Molec. Struct. (Theochem) , 179 (2006).[88] M. E. Casida, F. Gutierrez, J. Guan, F. Gadea, D. R. Salahub, and J. Daudey,
Charge-transfer correction for improved time-dependent local density approximationexcited-state potential energy curves: Analysis within the two-level model withillustration for H and LiH , J. Chem. Phys. , 7062 (2000).[89] M. E. Casida, Review: Time-dependent density-functional theory for moleculesand molecular solids , J. Mol. Struct. (Theochem) , 3 (2009).[90] M. E. Casida and M. Huix-Rotllant,
Progress in time-dependent density-functional theory , Annu. Rev. Phys. Chem. , in press.[91] M. E. Casida and D. R. Salahub, Asymptotic correction approach to improv-ing approximate exchange-correlation potentials: Time-dependent density-functional theory calculations of molecular excitation spectra , J. Chem. Phys., , 8918 (2000).[92] M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub,
Molecular excitationenergies to high-lying bound states from time-dependent density-functionalresponse theory: Characterization and correction of the time-dependent localdensity approximation ionization threshold , J. Chem. Phys. , 4439 (1998).[93] W. F. Chan, G. Cooper, and C. E. Brion,
Absolute optical oscillator strengths fordiscrete and continuum photoabsorption of carbon monoxide (7-200 ev) andtransition moments for the X Σ + → A Π system , Chem. Phys. , 123 (1993).[94] T. A. Arias, Multiresolution analysis of electronic structure: semicardinal andwavelet bases , Rev. Mod. Phys. , 267 (1999).[95] P. Fischer and M. Defranceschi, Looking at atomic orbitals through Fourier andwavelet transforms , Int. J. Quant. Chem. , 619 (1993).[96] J.-L. Calais, Wavelets–Something for quantum chemistry? , Int. J. Quant. Chem.58