Linear Scaling Density Matrix Real Time TDDFT: Propagator Unitarity \& Matrix Truncation
LLinear Scaling Density Matrix Real Time TDDFT: Propagator Unitarity & MatrixTruncation
Conn O’Rourke and David R. Bowler London Centre for Nanotechnology, University College London, 17-19 Gordon St,London, WC1H 0AH, United Kingdom a) (Dated: 5 November 2018) Real time, density matrix based, time dependent density functional theory proceedsthrough the propagation of the density matrix, as opposed to the Kohn-Sham orbitals.It is possible to reduce the computational workload by imposing spatial cut-off radiion sparse matrices, and the propagation of the density matrix in this manner providesdirect access to the optical response of very large systems, which would be otherwiseimpractical to obtain using the standard formulations of TDDFT. Following a briefsummary of our implementation, along with several benchmark tests illustrating thevalidity of the method, we present an exploration of the factors affecting the accuracyof the approach. In particular we investigate the effect of basis set size and matrixtruncation, the key approximation used in achieving linear scaling, on the propagatorunitarity and optical spectra. Finally we illustrate that, with an appropriate densitymatrix truncation range applied, the computational load scales linearly with thesystem size and discuss the limitations of the approach.Keywords: Density Matrix, TDDFT a) UCL Satellite, MANA International Centre for Materials Nanoarchitectonics (MANA), National Institutefor Materials Science (NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan; Electronic mail: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] F e b . INTRODUCTION Linear scaling or O ( N ) density functional theory (DFT), in which the computationalworkload scales linearly with the number of atoms in the system N , is now well established .In the standard approach to DFT, diagonalisation of an eigenvalue equation, or alternativelythe orthogonalisation of the Kohn-Sham states during minimisation of the energy, resultsin a severe computational bottleneck that limits the size of systems which can be studied.Working with the density matrix, upon which a truncation radius is applied, allows thecomputational workload to be made to scale linearly with N . Circumventing the size limi-tations of the standard approach in this manner allows vastly larger systems to be studied:for example calculations have now been performed on millions of atoms , in comparison tothe upper limit of around a thousand for the standard approach.While density functional theory is a ubiquitous tool in the arsenal of the electronic struc-ture theorist, it is limited to the study of ground-state properties. Extending DFT to thetime domain results in its excited state couterpart, time dependent density functional theory(TDDFT). Linear response TDDFT (LR-TDDFT), as developed by Casida , again suffersfrom a computational bottleneck which forces it to scale poorly with system size. LR-TDDFT requires the solution of an eigenvalue equation for a matrix written in the space ofelectron-hole pairs, which ostensibly scales as poorly as O ( N ). In practice this scaling canbe reduced, through efficient implementation and methods employing the Liouville-Lanczosapproach, to be as low as O ( N ) . For small systems LR-TDDFT is computationally fea-sible, and has been widely used, while for larger systems the scaling renders it unsuitable.It is also worth noting that linear scaling density matrix based LR-TDDFT, avoiding thepropagation of the density matrix, has also been recently demonstrated .An alternative approach to LR-TDDFT is the real time propagation of the time-dependent Kohn-Sham equations, pioneered by Yabana and Bertsch . Real time TDDFT(RT-TDDFT) proceeds by the construction of an effective Hamiltonian, followed by thedirect propagation of the Kohn-Sham orbitals using this Hamiltonian. Assuming both thenumber of occupied states (N KS ) and the number of mesh points (N M ) scale linearly withsystem size, RT-TDDFT will scale with the number of atoms, N, as N KS N M ∼ N . Asignificant prefactor in the form of the number of time steps and the computational effortfor construction of the Hamiltonian exists, making this method unsuitable for systems of2mall size. However the O (N ) scaling have made it the natural choice for tackling systemsof large size, and a complementary partner to Casida’s approach.In a similar manner to O ( N ) DFT, it is possible to improve upon the scaling of RT-TDDFT by propagating the density matrix, as opposed to propagating the Kohn-Shamorbitals directly. By applying a spatial truncation radius upon the density matrix, thecomputational workload can be reduced, opening up the possibility of studying excitedstates in large systems that cannot feasibly be examined with other methods. Although notwidely employed, this approach has been demonstrated to scale linearly with system size,and has been used to study several large systems; fullerene, sodium clusters, polyacetyleeoligomers, carbon nanotubes and silicon clusters to name a few .Several factors must be taken into consideration when employing this method, for exam-ple the accuracy of results produced will depend strongly on the range of truncation of thedensity matrix. Also when working in a non-orthogonal basis, as is the case in the CON-QUEST code, the overlap matrix will be well-ranged. However the inverse overlap, whichfeatures in the density matrix propagators, will not necessarily be. In order to ensure theunitarity of the propagation the propagtors must be carefully tested for matrix truncationerrors, and little discussion on the effect of matrix truncation upon propagator unitarityhave been presented elsewhere.In this paper we briefly summarize our implementation of RT-TDDFT in the CON-QUEST code, for completeness, and confirm its reliability. We then present several testsprobing the limitations of the method, and factors affecting accuracy. In particular we exam-ine the effect of matrix truncation, the key approximation used in achieving linear scaling,on the unitarity of the propagators used and optical spectra generated. II. COMPUTATIONAL APPROACH
Linear scaling approaches for excited state properties have existed for well over a decade(for a review see ), with the first approach utilising the locality inherent in the densitymatrix and being carried out at the semi-empirical level . Subsequent efforts again all tendto employ the nearsightedness of the density matrix, with the first full linear scaling TDDFTbeing done by Yam et. al. . Our approach follows that of Yam et. al. closely, with a fewdifferences; most notably we choose not to perform the orthogonalisation procedure via the3holesky decomposition and rather work in our non-orthogonal basis. As mentioned, linearscaling approaches to calculating the excited state properties in the frequency domain havealso been presented, by Yokojima et. al. , and more recently by Zuehlsdorff et. al. in theONETEP code . It is also worth noting that an approach for calculating the unoccupiedKohn-Sham states, via a basis optimisation approach which is also linear scaling, has alsobeen implemented in the ONETEP code .In a similar vein to the standard approaches to TDDFT in the time and frequency domain,the reformulations using the density matrix can be viewed as complementary to one another.The frequency domain approach is suitable for calculating the lowest optical excitations inthe system, but if the density matrix response involves higher excitations it will not besuitable. While the real-time density matrix approach employed here and by Yam et. al.calculates the full optical spectrum, it has a significant prefactor in the form of the numberof time steps needed for the numerical integration.In this section we briefly give an overview of the approach in our non-orthogonal basisset, and in the subsequent section we illustrate the effect of the basis set on the results, andthe reliability of the method with several tests on small molecules. A. Density Matrix RT-TDDFT
Rather than working with the conventional single particle Kohn-Sham orbitals, CON-QUEST works with the density matrix written in a seperable form in terms of a localisedbasis of support functions φ iα ρ ( r , r (cid:48) ) = (cid:88) iα,jβ φ iα ( r ) K iα,jβ φ jβ ( r (cid:48) ) (1)where φ iα is the α th support function centred on atom i . Support functions are a non-orthogonal basis set of localised orbitals, and have an overlap matrix given by: S α,β = (cid:90) φ iα ( r ) φ jβ ( r ) d r (2)Linear scaling behaviour can be obtained through applying a spatial cut-off on the densitymatrix. Beyond this cut-off radius the matrix elements are set to zero which, along withthe spatial limitation of the support functions, ensures that the number of non-zero density4atrix elements increases linearly with system size (for a fuller overview of the CONQUESTcode see ).RT-TDDFT is now well established , and implementations of density matrix RT-TDDFThave been reported elsewhere . Rather than employing an orthogonalisation procedurevia a Cholesky or L¨owdin decomposition, which will increase the range of the sparse matricesand is done elsewhere, we work in our non-orthogonal basis. Expanding the time-dependentKohn-Sham equations in this basis of non-orthogonal support functions, in the instancewhere the support functions are stationary with time, gives: i ∂∂t c ( t ) = S − Hc ( t ) (3)and i ∂∂t c † ( t ) = − c † ( t ) HS − (4)which describe the time dependence of the coefficients of our basis set expansion, c ( t ).This allows us to write the quantum Liouville equation of motion for our auxiliary densitymatrix K in the non-orthogonal support function basis: i ˙ K = S − HK − KHS − (5)The formal solution to this equation can be expressed as: K ( t ) = U ( t , t ) K ( t ) U † ( t , t ) (6)where U ( t,t ) is a propagator satisfying both: c ( t ) = U ( t , t ) c ( t ) (7) i ∂∂ t U ( t , t ) = S − H U ( t , t ) (8)Expressing the propagator in integral form we have: U ( t , t ) = T exp (cid:26) − i (cid:90) tt d τ S − H ( τ ) (cid:27) (9)where T is the time ordering operator. Evolution of the system for a total time, T = n ∆ t ,may be carried out piecewise in smaller intervals, allowing us to express the total evolutionoperator as the product of small time operators: U ( t, t ) (cid:39) N − (cid:89) n =0 U (( n + 1) ∆ t , n ∆ t ) (10)5here U ( t + ∆ t , t ) = exp (cid:2) − i S − H ( τ )∆ t (cid:3) (11)Evolution of the time dependent system is then reduced to the problem of approximatingthe propagator U ( t + ∆ t , t ). Two approximations exist in the definition of U ( t + ∆ t , t ),firstly that of approximating the matrix exponential exp( A ) and secondly the exact formof the matrix for which we wish to calulate the exponential. There are several methods forcalulating the exponential of a matrix , here we use the simplest approximation, a Taylorexpansion: exp( A ∆ t ) = I + ∞ (cid:88) n =1 ( A ∆ t ) n n ! (12)Similarly there are many different approaches for deciding which matrix exponential to useas a propagator. Three approximations have been implemented: the so called exponential-midpoint propagator (EM), the enforced time-reversal symmetry (ETRS) propagator andthe fourth order Magnus (M4) propagators, all of which are taken from the work of Marqueset al. on RT-TDDFT propagators, and are briefly described in our non-orthogonal basisfor completeness.The exponential midpoint propagator approximates the U ( t + ∆ t, t ) by the exponentialtaken at τ = t + ∆ t/ U EM ( t + ∆ t , t ) = exp (cid:26) − i S − H (cid:18) t + ∆ t (cid:19)(cid:27) (13)Implicitly enforcing time-reversibility, such that propagating forward from t and back-wards from t + ∆ t by ∆ t / enforced time-reversal symmetry method: U ET RS ( t + ∆ t , t ) = exp (cid:26) − i ∆ t S − H ( t + ∆ t ) (cid:27) × exp (cid:26) − i ∆ t S − H ( t ) (cid:27) (14)Using the Magnus operator the exponential solution to Schr¨odinger equation for a time-dependent Hamiltonian may be written as : U M ( t + ∆ t , t ) = exp { M G } (15)6here M G4 is an infinite series of integrals providing an exact solution. Truncating thisexpansion to fourth order and approximating the integrals using Gauss-Legendre points asin gives in our non-orthogonal basis: M G = − i ∆ t (cid:2) S − H ( t ) + S − H ( t ) (cid:3) − √ t (cid:2) S − H ( t ) , S − H ( t ) (cid:3) (16)where t , = t + [1 / ± √ / t .It is important to note the presence of the inverse overlap matrix S − in these propaga-tors, and again consider that while the overlap matrix will be well-ranged and suitable fortruncation, the inverse overlap is not necessarily so. We therefore need to carefully test thesparsity of the product S − H , and its effect on the unitarity of our propagators. B. Linear Response
The idea behind extracting optical transitions from the linear response of a system to anexternal electric field is well known . Propagating in real time provides direct access tothe time-dependent charge density, and therefore the electronic response to external fields.Applying a time dependent external electric field polarised along axis j , δ v ext ( r , t ) = − E j ( t ) · r allows us to examine the time-dependent response of the system. Application of this electricfield will produce an induced time-dependent dipole moment: P ( t ) = P (0) − (cid:90) d r n ( r , t ) r . (17)As an example of the calculated repsonse of a system to an applied electric field, figure 1illustrates the induced dipole response of a benzene molecule on application of a field witha Gaussian time profile, centred at t = 0.Access to the time-dependent dipole moment allows us to calculate the time dependentpolarisability: α ij ( ω ) = (cid:82) dt e iωt P i ( t ) (cid:82) dt e iω t E j ( t )The imaginary part of the polarisability is directly proportional to the absorption crosssection, σ ( ω ) and the experimentally observed strength function, S ( ω ).7 IG. 1. Applied electric field and induced dipole moment for a benzene molecule. (∆ t = 0 .
03 a.u. ≈ S ( ω ) = 2 ωπ Im (cid:18) T r ( α µj ( ω )) (cid:19) (18)As noted by Tsolakidis et. al., the approach satisifies the f-sum rule and the integrationof the strength function over energy gives the number of electrons, which may be used as ameasure of the completeness of the basis set .Density matrix RT-TDDFT therefore has the potential to be an extremely useful tool fortheoretically predicting the electronic absorption spectra of large system. III. SMALL MOLECULES
In order to verify that our implementation is correct we have performed tests on severalsystems for which the electronic transitions have been studied experimentally and theoret-ically elsewhere, allowing us to make direct comparisons. For this purpose we have chosenfour small molecules (Carbon monoxide, Methane, Ethylene and Benzene) and used our im-plementation to calculate the optical absorption spectra within the TDLDA approximation.Meaningful comparison of our results with experiment requires the identification the elec-tronic transitions to which the peaks in our calculated absorption spectra correspond. As wehave mentioned in Casida’s approach information about electronic transitions is inherentlyproduced, while in RT-TDDFT it is not. 8 ransition Basis Set Ref Expt. π → π ∗ π → s π → π ∗ π → s π → π ∗ ) and Rydberg ( π → s ) excitations for the C H molecule. It is often possible to identify the corresponding transition by examining the polarisationand energy of peaks and comparing to that of optically allowed transitions experimentally.Where possible, in order to more confidently assign peaks of our calculated absorptionspectra to particular electronic transitions, we have followed the procedure in wherebya sinusoidal electric field tuned to a particular excitation mode is applied. A resultingelectronic resonance is set up, allowing us to examine the difference between ground statecharge density and excited state charge density and thereby infer the electronic transition. A. Basis Sets
Our support functions are expanded in a basis of numerical orbitals, in this case pseudo-atomic orbitals generated following the approach of the Siesta code . These PAOs areeigenfunctions of the atomic pseudopotentials with a confinement energy shift used to de-termine a radial cut-off for the orbitals, beyond which they are zero. This confinementenergy provides a single parameter to define the cut off radii for different orbitals, and is theenergy each orbital obtains on being confined by an infinite potential to a particular radius.It is clear that a minimal basis with which ground state properties are accurately reproducedwill generally not be satisfactory for calculating excited state properties, and therefore weillustrate the basis set dependence of two selected transitions for the C H molecule.Multiple orbitals per angular momentum channel can be used (multiple- ζ ), with the shape9 olecule Transition RT-TDDFT Ref Expt(eV.)CO σ → π ∗ CH T →
3s 9.22 9.27 C H π → π ∗ C H π → π ∗ ∼ TABLE II. Comparison of calculated TDLDA transition energies for small molecules with othervalues and experiment. Conquest results obtained with 5Z4P basis sets, with the exception ofbenzene (2Z2P). of multiple orbitals determined by a split norm procedure . This procedure uses a parameterto define the norm of a numerical orbital outside some radius where they match the tail ofthe first zeta PAO, and within this radius the vary smoothly to the origin. Subtractingthis numerical orbital from the original PAO gives the multiple-zeta orbital. Of course it ispossible to define these radii by hand and fine tune the basis set. In addition to multiplezeta, polarisation orbitals can be included within the basis set, and are obtained by solvingthe same pseudo-atomic problem but with an applied electric field.We use the notation SZ, 2Z, 3Z, 4Z to describe single zeta, double zeta, triple zeta and soon. Similarly we describe the number of polarisation orbitals included in the basis by SZP,SZ2P and SZ3P (one, two and three polarised orbitals respectively).To first gauge the effect of varying our basis set on the results we have performed calcula-tions on the ethylene molecule with varying numbers of PAOs and two different confinementenergies. The basis sets have been generated with a confinement energy of 1 meV and 5meV, resulting in confinement radii of 4.93 and 4.24 ˚A for the carbon atoms respectively,and radii of 4.77 and 4.21 ˚A for the hydrogen atoms respectively. The total run time was14.51 fs. (600 au.) with a time step of ∼ π → s transition show a wide variation with basis set choice,while the π → π ∗ valence transition varies less. This is in line with expectation, given themore diffuse nature of the Rydberg transition we would expect its description to requirelarger basis. The effect of systematically increasing the number of basis functions is to10mprove our results with respect to that of the reference values. Similarly increasing thecut-off radii, by reducing the confinement energy, tends to improve the quality of the result.This is to be expected, as increasing the size of our basis set, while systematically increasingthe range, will maximise the variational degrees of freedom available to describe our timedependent density matrix.However our values are still far from those computed elsewhere, and we find generallythat for small molecules it is essential to use a large basis with multiple extended polarisationorbitals in order to produce results in line with other works. In addition we find that finetuning the radial cut-offs by hand, as opposed to using the confinement energy and splitnorm procedure, can allow us to improve the quality of our results for small molecules. B. Small Molecule Results
Exhibited in table II are the calculated transitions for our four test molecules. In thecase of the smallest molecules (carbon monoxide, ethylene, and methane) a hand tuned5Z4P basis set is employed, while for benzene the result is obtained using a 2Z2P basiswith a 5meV confinement energy (all the calculations satisfy the f-sum rule to > IV. PROPAGATOR UNITARITY
Having demonstrated the correctness of our implementation and explored the influenceof basis sets, we now turn to our main concern, the effects of localisation in linear scaling11 i)(ii)
FIG. 2. (i) : Absorption strength function for carbon monoxide from RT-TDDFT and experiment.Experimental data taken from . (ii) Absorption strength function for Benzene from RT-TDDFT.Experimental data taken from . methods on the accuracy of results.We wish the total charge in our system to remain stable, and in order for this to be thecase the propagators must be unitary with respect to the non-orthogonal basis set: U † U − I = 0 (19)where U is our propagator matrix and I is the identity matrix.From our approximation for the matrix exponential, eq. 12, it can be shown that, if itwere exact, our propagators would indeed exhibit this property. However, as it is impossiblefor us to store an infinite sum on our computer, we must truncate our Taylor expansionat some point. Doing so will introduce errors, with two factors affecting the scale of the12 IG. 3. Plot of the absolute values of matrix U † U − I (on a base 10 log scale), illustrating thepropagator unitarity for the exponential midpoint propagator, for varying time step sizes (in a.u.).The system studied is a single benzene molecule, and the matrix is shown at the end of a 10 a.u.run. break from unitarity; the time step and the number of terms in our summation. Whilewe can extend our expansion arbitrarily, and reduce the time step arbitrarily, we wish toavoid excess computational expense by keeping the expansion as small as possible and thetime step as large as possible within some acceptable margin of accuracy. We can directlyexamine the unitarity of our propagators through equation 19. A. Time-Step Dependence
As a test we have examined the extent of the break from unitarity for a range of time-stepsand number of terms in the matrix exponential expansion. We have used a small moleculefor the purpose, benzene, with a small applied electric field perturbation with a Gaussianprofile centered on t = 0.Exhibited in figure 3 we can see the dependence on simulation time step of the propaga-tor unitarity, with the obvious trend being that as the time step is reduced the propagatorapproaches unitarity. We can see that even for time steps up to ∼ .
15 a.u. the propagatormaintains its unitarity to a high degree (similar results were obtained for each of the prop-agators). The corresponding effect on the charge conservation can be seen in figure 4 and,as expected we see that as the time step increases the conservation of charge deteriorates13ith the propagation eventually becoming unstable for large timesteps. While the maxi-mum permissible timestep will depend on the system under study, we found that generallya timestep of 0.06 a.u. or below provided satisfactory charge conservation.The form of our propagators requires the extrapolation of the Hamiltonian matrix tosome unknown point beyond the current time t, H + . As suggested by Marques et al. inorder to minimise errors it is possible to carry this procedure out self-consistently. In ourcase meaning that we propagate K ( t ) to K ( t + ∆ t ) based on an extrapolated Hamiltonian.We then construct a new Hamiltonian matrix H ( t + ∆ t ) using K ( t + ∆ t ). H + can then beinterpolated from Hamiltonian matrices for times up to and including ( t +∆ t ), and the wholeprocedure is iterated until some self-consistency criteria is obtained. Generally speaking thisprocedure is performed three times in the early stages of a run, following a perturbation,and reduces to two as the run progresses. The effect of not performing this self-consistencyprocedure on the charge conservation can be seen in figure 4. While the self-consistency cycleis found to improve the charge conservation, in reality for small time steps the differencein charge conservation and calculated properties is not found to be significant enough towarrant the extra computational load of constructing the Hamiltonian matrix several timesper time-step. As a compromise we enforce the self-consistency only for a small number ofsteps ( ∼ − U M4 , is advantageous , butfor our present work this is not the case and we have opted for the simplest exponentialmidpoint propagator throughout. B. Matrix Exponential Truncation
The effect of truncating the Taylor expansion used to evaluate the matrix exponential onthe unitarity of the propagator can be seen in figure 6. We see that reducing the number ofterms reduces the unitarity of the propagator, as expected. Looking at figure 5 the conver-gence of the charge conservation with the number of terms in the exponential expansion can14
IG. 4. Variation in total charge (on a base 10 log scale) with time step size, following a 10 au. runfor benzene using all three propagators. Also included is charge variation for the EM propagatorwithout the self-consistent propagator step (see text for details).!FIG. 5. Absolute variation in total charge (on a base 10 log scale) with the number of terms in ourmatrix exponential expansion, following a 20 au. run for benzene using the EM propagator with atime step of 0.04 au. be seen. We find that we reach good convergence with six terms included in the expansion,and we opt for this level of accuracy throughout the remainder of the paper.15
IG. 6. Plot of the absolute values of the matrix U † U − I (on a base 10 log scale), illustratingthe propagator unitarity for the exponential midpoint propagator, for differing number of termsin the Taylor expansion for our propagator. The system studied is a single benzene molecule, andthe matrix is shown at the end of a 10 au. run (dt = 0.04 au.) V. ALKANE MOLECULES: TESTING MATRIX TRUNCATION EFFECTS
In this section we perform calculations on long chain alkane molecules.Our aim is to exam-ine the effect of matrix truncation on the propagation of the density matrix and propagatorunitarity, along with the computational scaling with system size.As a first step we calculate the absorption spectra for the C H molecule for severaldifferent basis sets using the generalised gradient PBE functional (all further calculationsin this section are performed with this functional), and the results can be seen in figure 7.Experimentally as the length of the alkane carbon chain increases, the absorption onset isfound to reduce, and the reported adsorption onset for C H is ∼
175 nm. ( ∼ ∼ , we would expect the addition of more diffusePAOs to improve the description of these excitations. Given the well documented difficulties16 i) (ii) FIG. 7. Basis set variation of the calculated alkane optical absorption spectra. (i)
Effect of in-creasing the number of PAOs in the basis set and (ii) the effect of extending the radii of the basisfunctions are shown. of TDDFT to accurately describe Rydberg transitions , and given that this is not our aimin any case, we proceed to carry out our tests with the SZP and SZ2P basis sets generatedusing a confinement energy of 55meV (radial cut off for the PAOs is 3.31˚A and 3.12˚A forcarbon and hydrogen respectively).Yam et al. have previously studied the long chain alkanes within the linear scalingexcited state regime , calculating the absorption onset at around 8 eV for C H with theLDA functional. However little discussion of the effects of matrix truncation on propagatorunitarity have been presented elsewhere. A. Propagator Truncation
The use of a basis of non-orthogonal atomic orbitals requires the inverse overlap matrixfor our propagation (indeed this matrix is required for ground state calculations in anycase), as seen in equation 11. In order to compute the inverse overlap matrix Conquestuses Hotelling’s method , however for poorly conditioned overlap matrices computing theinverse overlap matrix can prove difficult. In our current implementation of TDDFT theatoms remain stationary and so too, therefore, does the overlap matrix. Therefore we17 IG. 8. Average absolute error in the S -1 H (left) and S -1 (right) matrix elements with matrix rangefor the C H molecule. SZP basis set is used, generated with a 55meV confinement potential. have also included the possibility of computing the inverse overlap with the SCALAPACKroutines. Although the scaling will not be linear, computing the inverse overlap in this waymakes only a relatively small contribution to our total TDDFT runtime, as we only calculatethe inverse overlap once at t = 0.While it is apparent that the overlap matrix will be sparse, allowing it to be truncated,the inverse of a sparse matrix will not in general be sparse itself. We have therefore testedthe effect of truncating both the S -1 matrix and the S -1 H matrix on the propagation. Figure8 shows the average absolute error in the matrix elements of S -1 and the S -1 H matricescaused by truncation (the error in S -1 elements given is the average of the elements of the S -1 S - I matrix, and the error in the S -1 H is calculated with the values from an untruncated S -1 matrix).As the range of the matrices increases the error caused by the truncation convergestowards zero, as we expect. The S -1 matrix converges less quickly than the S -1 H matrix,indicating that it is more dense than the S -1 H matrix. The effect truncation of thesematrices has on the unitarity of the propagators can be seen in figure 9. We see that theunitarity converges as the S -1 H range increases, and the propagators are converged with arange of around ∼ S -1 H matrix is indeed sparse,while the S -1 matrix is less so, and we can safely truncate it. It is important to not thatwe don’t explicitly use the S -1 in our propagators, only the S -1 H matrix. Although it18 IG. 9. Plot of the absolute values of the matrix U † U − I (on a base 10 log scale), illustrating thepropagator unitarity for differing truncation ranges of the S -1 and S -1 H matrices for the C H molecule makes sense to truncate the S -1 matrix, given that we are truncating S -1 H and that theHamiltonian matrix is sparse. We can see this by noting that the unitarity of the propagatorin figure 9 is also well converged for each of the truncation ranges imposed on the inverseoverlap.As additional atoms are added the Hamiltonian matrix, overlap matrix, and the inverseoverlap will vary. Increasing the system size may therefore affect the ranges of these matrices.While we only use the S -1 H matrix in our calculation, comparison of the density of bothmatrices have been included. We have tested this effect by fixing the S -1 and S -1 H rangesat 30 and 35 Bohr respectively, and examined the error in the truncated S -1 H matrix withsystem size with the results shown in figure 10. We see that the error changes slightlyon increasing system size, but converges as the size increases. Consequently the propagatorunitarity was found to exhibit the same trend. This illustrates that the S -1 H is well ranged,irrespective of system size, allowing us to impose a cut-off radii on both of these matrices.In effect this ensures that as the system size increases, the computational load can be madeto scale linearly.Similarly, increasing the number of basis functions will directly affect the overlap matrix,and consequently the inverse overlap and the propagator. In order to gauge the extent ofthis effect we have examined the C H molecule with a larger basis set (SZ2P as opposed19 IG. 10. Average absolute error in the S -1 H matrix elements with system size.FIG. 11. Average value of the U † U − I matrix with S -1 H matrix range for the C H molculecalculated with a SZ2P basis set. to SZP). Exhibited in figure 11 is the absolute value of the U † U − I matrix with S -1 H matrix truncation range. Despite the larger number of basis set functions we see that the S -1 H matrix is still well ranged, although the range is wider when compared to the SZPresults of figure 9, and again a truncation will lead to a computational load that scaleslinearly with system size.A further point to note is that it is possible to avoid the use of the inverse overlapmatrix in the TDDFT propagation altogether. Yam et al. have employed a Choleskyorthogonalisation scheme to bypass the need for the inverse overlap .n However using thisscheme requires the inverse of the Cholesky decomposition, and it is not apparent that it will20 IG. 12. K matrix truncation radii dependence: Spectra generated for the C H molecule atvarying density matrix cut-off radii. (Total run time of 400 a.u. at a time step of 0.05 a.u.) be more sparse than the inverse overlap. It is possible that this scheme might improve thecalculation of the propagator, as the orthogonalised Hamiltonian may be more localised thanour S -1 H matrix. Calculating the Cholesky decomposition can be made to scale linearly,and implementation of this alternative method has already begun in order to contrast thetwo approaches. However the parallelisation of Cholesky inversion is difficult given theConquest matrix storage, and inversion of the overlap matrix remains important. B. Density Matrix Truncation, Scaling and Limits
Finally we examine the effect of truncating the density matrix, and have performedcalculations generating spectra for the C H molecule at varying truncation radii, R Cut ,of the density matrix. Typically for ground state calculations a suitable typical densitymatrix truncation range is around 16-20 Bohr. The results can be seen in figure 12, andgenerally we find that as the density matrix cut-off increases the spectra tend to converge,as expected, with higher lying states requiring a larger cut-off to converge. We can see fromthe comparison of R Cut = 30 and R Cut = 35 that there is good agreement for the initialtransitions, as well as the general shape of the spectra.21pplying this R Cut = 35 Bohr cut-off (along with a cut off of 35 Bohr. on the S − H ma-trix) we can examine the computational scaling with system size, with the results exhibitedin figure 13. Clear linear scaling of the computational workload up to well over 1000 atomsis exhibited, illustrating the potential power of the method.Finally a few comments on the limits of the approach must be made. TDDFT forlong-range charge transfer is well known to be poorly described by local and semi-localfunctionals . While we have employed LDA and GGA functionals here, linear scaling ex-act exchange has also been recently implemented in the Conquest code, allowing the use ofnon-local functionals with this approach in the future.While the near-sightedness principle dictates that the ground-state density matrix isexponentially localised for well gapped systems, there is no formal justification for the lo-calisation of the response density matrix. As noted in , for systems with well localisedexcitations it would be expected that the response density matrix could be truncated safelyand linear scaling achieved, while for systems with delocalised excitations this will not bethe case. FIG. 13. Computational TDDFT run time versus system size for long chain alkane molecules. Thesystem was run with a timestep of 0.05 a.u. for a total time of 10 a.u. A matrix truncation range, R Cut = 35 a.u., has been applied. I. CONCLUSIONS
We have outlined our implementation of real-time time dependent density functional the-ory in the Conquest O (N) code. We have demonstrated the soundness of the implementationthrough benchmark tests for small molecules, and also discussed the effect of basis set andsystem sizes on the results. O ( N ) approaches utilise the density matrix, as opposed to working directly with Kohn-Sham orbitals, providing a route to the linear scaling computational time with system sizeby its truncation. We have discussed the range of our propagator matrices for an alkanechain test system, and the implications of this matrix truncation on the unitarity of thepropagation. Similarly we have examined the effect of truncating the density matrix onthe calculated optical absorption spectra, showing that the range required is much moreextended than that required for converged ground state properties. Nevertheless, we haveshown that accurate linear scaling TDDFT calculations are practical. While the impactof localisation cut-off in the charge density matrix on these TDDFT calculations is a topicwarranting further study, we have shown that in truncating these matrices at a suitablepoint we obtain a computational load that increases linearly with system size. This offersa complementary approach to the usual Casida linear response approach: linear responseTDDFT is well suited to relatively small systems, while linear scaling RT-TDDFT offersa viable method for studying excitations in large systems. We have shown linear scalingbeyond 1,000 atoms, and 10,000+ atoms are perfectly practical with the excellent parallelscaling available in Conquest. ACKNOWLEDGEMENTS
C.O’R. is supported by the MANA-WPI project and D.R.B. was funded by the RoyalSociety. We thank Umberto Terranova for useful discussions. This work made use of thefacilities of ARCHER, the UK’s national high-performance computing service, which isprovided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, andfunded by the Office of Science and Technology through EPSRC’s High End ComputingProgramme. Calculations were performed at HECToR through the UKCP Consortium.The authors acknowledge the use of the UCL Legion High Performance Computing Facility,23nd associated support services, in the completion of this work.
REFERENCES D. R. Bowler and T. Miyazaki, “ O (n) methods in electronic structure calculations,” Re-ports on Progress in Physics , 036503 (2012). D. R. Bowler and T. Miyazaki, “Calculations for millions of atoms with density functionaltheory: linear scaling shows its potential,” Journal of Physics: Condensed Matter ,074207 (2010). J. VandeVondele, U. Bortnik, and J. Hutter, “Linear scaling self-consistent field calcu-lations with millions of atoms in the condensed phase,” Journal of Chemical Theory andComputation , 3565–3573 (2012), http://dx.doi.org/10.1021/ct200897x. M. E. Casida, in Recent Developments and Applications of Modern Density Functional Theory(Elsevier, Amsterdam, 1996). L. Kronik, A. Makmal, M. L. Tiago, M. M. G. Alemany, M. Jain, X. Huang, Y. Saad,and J. R. Chelikowsky, “Parsec the pseudopotential algorithm for real-space electronicstructure calculations: recent advances and novel applications to nano-structures,” physicastatus solidi (b) , 1063–1079 (2006). D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, “Turbo charging time-dependent density-functional theory with lanczos chains,” The Journal of Chemical Physics , 154105(2008). T. J. Zuehlsdorff, N. D. M. Hine, J. S. Spencer, N. M. Harrison, D. J. Riley, and P. D.Haynes, “Linear-scaling time-dependent density-functional theory in the linear responseformalism,” The Journal of Chemical Physics , 064104 (2013). K. Yabana and G. F. Bertsch, “Time-dependent local-density approximation in real time,”Phys. Rev. B , 4484–4487 (1996). S. Yokojima and G. Chen, “Time domain localized-density-matrix method,” ChemicalPhysics Letters , 379 – 383 (1998). S. Yokojima and G. Chen, “Linear scaling calculation of excited-state properties of poly-acetylene,” Phys. Rev. B , 7259–7262 (1999). S. Yokojima and G. Chen, “Linear-scaling localized-density-matrix method for the groundand excited states of one-dimensional molecular systems,” Chemical Physics Letters ,2440 – 544 (1999). S. Yokojima, D. Zhou, and G. Chen, “Linear-scaling computation of ground state withtime-domain localized-density-matrix method,” Chemical Physics Letters , 495 – 498(1999). W. Liang, S. Yokojima, and G. Chen, “Generalized linear-scaling localized-density-matrixmethod,” The Journal of Chemical Physics , 1844–1855 (1999). W. Liang, S. Yokojima, D. Zhou, and G. Chen, “Localized-density-matrix method and itsapplication to carbon nanotubes,” The Journal of Physical Chemistry A , 2445–2453(2000), http://pubs.acs.org/doi/pdf/10.1021/jp990818a. A. Tsolakidis, D. S´anchez-Portal, and R. M. Martin, “Calculation of the optical responseof atomic clusters using time-dependent density functional theory and local orbitals,” Phys.Rev. B , 235416 (2002). C. Yam, Q. Zhang, F. Wang, and G. Chen, “Linear-scaling quantum mechanical methodsfor excited states,” Chem. Soc. Rev. , 3821–3838 (2012). C. Y. Yam, S. Yokojima, and G. Chen, “Localized-density-matrix implementation of time-dependent density-functional theory,” The Journal of Chemical Physics , 8794–8803(2003). C. Yam, S. Yokojima, and G. Chen, “Linear-scaling time-dependent density-functionaltheory,” Phys. Rev. B , 153105 (2003). L. E. Ratcliff, N. D. M. Hine, and P. D. Haynes, “Calculating optical absorption spectrafor large systems using linear-scaling density functional theory,” Phys. Rev. B , 165131(2011). C. Goringe, E. Hernndez, M. Gillan, and I. Bush, “Linear-scaling dft-pseudopotential cal-culations on parallel computers,” Computer Physics Communications , 1 – 16 (1997). G. Moore, “Orthogonal polynomial expansions for the matrix exponential,” Linear Algebraand its Applications , 537 – 559 (2011), ¡ce:title¿Special Issue: Dedication to PeteStewart on the occasion of his 70th birthday¡/ce:title¿. A. Castro, M. A. L. Marques, and A. Rubio, “Propagators for the time-dependent kohn–sham equations,” The Journal of Chemical Physics , 3425–3433 (2004). W. Magnus, “On the exponential solution of differential equations for a linear operator,”Communications on Pure and Applied Mathematics , 649–673 (1954).25 N. N. Matsuzawa, A. Ishitani, D. A. Dixon, and T. Uda, “Time-dependent den-sity functional theory calculations of photoabsorption spectra in the vacuum ul-traviolet region,” The Journal of Physical Chemistry A , 4953–4962 (2001),http://pubs.acs.org/doi/pdf/10.1021/jp003937v. V. J. Hammond and W. C. Price, “Oscillator strengths of the vacuum ultra-violet absorp-tion bands of benzene and ethylene,” Trans. Faraday Soc. , 605–610 (1955). S. K. Min, Y. Cho, and K. S. Kim, “Efficient electron dynamics with the planewave-basedreal-time time-dependent density functional theory: Absorption spectra, vibronic elec-tronic spectra, and coupled electron-nucleus dynamics,” The Journal of Chemical Physics , 244112 (2011). C. Hu, O. Sugino, and Y. Miyamoto, “Modified linear response for time-dependent density-functional theory: Application to rydberg and charge-transfer excitations,” Phys. Rev. A , 032508 (2006). W. Chan, G. Cooper, and C. Brion, “Absolute optical oscillator strengths for discreteand continuum photoabsorption of carbon monoxide (7-200 ev) and transition momentsfor the x σ + → a π system,” Chemical Physics , 123 – 138 (1993). N. N. Matsuzawa, A. Ishitani, D. A. Dixon, and T. Uda, “Time-dependent den-sity functional theory calculations of photoabsorption spectra in the vacuum ul-traviolet region,” The Journal of Physical Chemistry A , 4953–4962 (2001),http://pubs.acs.org/doi/pdf/10.1021/jp003937v. M. G. Curtis and I. C. Walker, “Low-energy electron-impact excitation of methane, silane,tetrafluoromethane and tetrafluorosilane,” J. Chem. Soc., Faraday Trans. 2 , 659–670(1989). K. Yabana and G. F. Bertsch, “Time-dependent local-density approximation in real time:Application to conjugated molecules,” International Journal of Quantum Chemistry ,55–66 (1999). E. Koch and A. Otto, “Optical absorption of benzene vapour for photon energies from 6ev to 35 ev,” Chemical Physics Letters , 476 – 480 (1972). J. M. Soler, E. Artacho, J. D. Gale, A. Garca, J. Junquera, P. Ordejn, and D. Snchez-Portal, “The siesta method for ab initio order- n materials simulation,” Journal of Physics:Condensed Matter , 2745 (2002). 26 J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation madesimple,” Phys. Rev. Lett. , 3865–3868 (1996). E. A. Costner, B. K. Long, C. Navar, S. Jockusch, X. Lei, P. Zimmerman, A. Campion,N. J. Turro, and C. G. Willson, “Fundamental optical properties of linear and cyclicalkanes: Vuv absorbance and index of refraction,” The Journal of Physical Chemistry A , 9337–9347 (2009), pMID: 19630422, http://pubs.acs.org/doi/pdf/10.1021/jp903435c. M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, “Molecular excitationenergies to high-lying bound states from time-dependent density-functional response the-ory: Characterization and correction of the time-dependent local density approximationionization threshold,” The Journal of Chemical Physics , 4439–4449 (1998). A. H. R. Palser and D. E. Manolopoulos, “Canonical purification of the density matrix inelectronic-structure theory,” Phys. Rev. B , 12704–12711 (1998). A. Dreuw, J. L. Weisman, and M. Head-Gordon, “Long-range charge-transfer excitedstates in time-dependent density functional theory require non-local exchange,” The Jour-nal of Chemical Physics119