# Linear Scaling Density Matrix Real Time TDDFT: Propagator Unitarity & Matrix Truncation

###### Abstract

Real time, density matrix based, time dependent density functional theory proceeds through 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 radii on sparse matrices, and the propagation of the density matrix in this manner provides direct access to the optical response of very large systems, which would be otherwise impractical to obtain using the standard formulations of TDDFT. Following a brief summary of our implementation, along with several benchmark tests illustrating the validity of the method, we present an exploration of the factors affecting the accuracy of the approach. In particular we investigate the effect of basis set size and matrix truncation, the key approximation used in achieving linear scaling, on the propagator unitarity and optical spectra. Finally we illustrate that, with an appropriate density matrix truncation range applied, the computational load scales linearly with the system size and discuss the limitations of the approach.

UCL Satellite, MANA] International Centre for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan

UCL Satellite, MANA] International Centre for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan

## I Introduction

Linear scaling or density functional theory (DFT), in which the computational workload scales linearly with the number of atoms in the system , is now well established Bowler and Miyazaki (2012). In the standard approach to DFT, diagonalisation of an eigenvalue equation, or alternatively the orthogonalisation of the Kohn-Sham states during minimisation of the energy, results in 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 the computational workload to be made to scale linearly with . Circumventing the size limitations 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 Bowler and Miyazaki (2010); VandeVondele, Bor¿tnik, and Hutter (2012), in comparison to the upper limit of around a thousand for the standard approach.

While density functional theory is a ubiquitous tool in the arsenal of the electronic structure theorist, it is limited to the study of ground-state properties. Extending DFT to the time domain results in its excited state couterpart, time dependent density functional theory (TDDFT). Linear response TDDFT (LR-TDDFT), as developed by CasidaCasida (1996), again suffers from 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 of electron-hole pairs, which ostensibly scales as poorly as . In practice this scaling can be reduced, through efficient implementation and methods employing the Liouville-Lanczos approach, to be as low as Kronik et al. (2006); Rocca et al. (2008). For small systems LR-TDDFT is computationally feasible, 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 the propagation of the density matrix, has also been recently demonstratedZuehlsdorff et al. (2013).

An alternative approach to LR-TDDFT is the real time propagation of the time-dependent Kohn-Sham equations, pioneered by Yabana and BertschYabana and Bertsch (1996). Real time TDDFT (RT-TDDFT) proceeds by the construction of an effective Hamiltonian, followed by the direct propagation of the Kohn-Sham orbitals using this Hamiltonian. Assuming both the number of occupied states (N) and the number of mesh points (N) scale linearly with system size, RT-TDDFT will scale with the number of atoms, N, as . A significant prefactor in the form of the number of time steps and the computational effort for construction of the Hamiltonian exists, making this method unsuitable for systems of small size. However the (N) scaling have made it the natural choice for tackling systems of large size, and a complementary partner to Casida’s approach.

In a similar manner to DFT, it is possible to improve upon the scaling of RT-TDDFT by propagating the density matrix, as opposed to propagating the Kohn-Sham orbitals directly. By applying a spatial truncation radius upon the density matrix, the computational workload can be reduced, opening up the possibility of studying excited states in large systems that cannot feasibly be examined with other methods. Although not widely employed, this approach has been demonstrated to scale linearly with system size, and has been used to study several large systems; fullerene, sodium clusters, polyacetylee oligomers, carbon nanotubes and silicon clusters to name a fewYokojima and Chen (1998, 1999a, 1999b); Yokojima, Zhou, and Chen (1999); Liang, Yokojima, and Chen (1999); Liang et al. (2000); Tsolakidis, Sánchez-Portal, and Martin (2002).

Several factors must be taken into consideration when employing this method, for example the accuracy of results produced will depend strongly on the range of truncation of the density matrix. Also when working in a non-orthogonal basis, as is the case in the CONQUEST code, the overlap matrix will be well-ranged. However the inverse overlap, which features in the density matrix propagators, will not necessarily be. In order to ensure the unitarity of the propagation the propagtors must be carefully tested for matrix truncation errors, and little discussion on the effect of matrix truncation upon propagator unitarity have been presented elsewhere.

In this paper we briefly summarize our implementation of RT-TDDFT in the CONQUEST code, for completeness, and confirm its reliability. We then present several tests probing the limitations of the method, and factors affecting accuracy. In particular we examine 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 Yam et al. (2012)), with the first approach utilising the locality inherent in the density matrix and being carried out at the semi-empirical levelYokojima and Chen (1998). Subsequent efforts again all tend to employ the nearsightedness of the density matrix, with the first full linear scaling TDDFT being done by Yam et. al. Yam, Yokojima, and Chen (2003a, b). Our approach follows that of Yam et. al. closely, with a few differences; most notably we choose not to perform the orthogonalisation procedure via the Cholesky decomposition and rather work in our non-orthogonal basis. As mentioned, linear scaling approaches to calculating the excited state properties in the frequency domain have also been presented, by Yokojima et. al.Yokojima and Chen (1998, 1999b), and more recently by Zuehlsdorff et. al. in the ONETEP code Zuehlsdorff et al. (2013). It is also worth noting that an approach for calculating the unoccupied Kohn-Sham states, via a basis optimisation approach which is also linear scaling, has also been implemented in the ONETEP codeRatcliff, Hine, and Haynes (2011).

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 in the system, but if the density matrix response involves higher excitations it will not be suitable. 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 number of time steps needed for the numerical integration.

In this section we briefly give an overview of the approach in our non-orthogonal basis set, and in the subsequent section we illustrate the effect of the basis set on the results, and the reliability of the method with several tests on small molecules.

### ii.1 Density Matrix RT-TDDFT

Rather than working with the conventional single particle Kohn-Sham orbitals, CONQUEST works with the density matrix written in a seperable form in terms of a localised basis of support functions

(1) |

where is the support function centred on atom i. Support functions are a non-orthogonal basis set of localised orbitals, and have an overlap matrix given by:

(2) |

Linear scaling behaviour can be obtained through applying a spatial cut-off on the density matrix. Beyond this cut-off radius the matrix elements are set to zero which, along with the spatial limitation of the support functions, ensures that the number of non-zero density matrix elements increases linearly with system size (for a fuller overview of the CONQUEST code see Goringe et al. (1997)).

RT-TDDFT is now well establishedYabana and Bertsch (1996), and implementations of density matrix RT-TDDFT have been reported elsewhereYam et al. (2012); Yam, Yokojima, and Chen (2003b). Rather than employing an orthogonalisation procedure via a Cholesky or Löwdin decomposition, which will increase the range of the sparse matrices and is done elsewhere, we work in our non-orthogonal basis. Expanding the time-dependent Kohn-Sham equations in this basis of non-orthogonal support functions, in the instance where the support functions are stationary with time, gives:

(3) | ||||

and | ||||

(4) |

which describe the time dependence of the coefficients of our basis set expansion, . This allows us to write the quantum Liouville equation of motion for our auxiliary density matrix K in the non-orthogonal support function basis:

(5) |

The formal solution to this equation can be expressed as:

(6) |

where is a propagator satisfying both:

c(t) | (7) | |||

(8) |

Expressing the propagator in integral form we have:

(9) |

where is the time ordering operator. Evolution of the system for a total time, , may be carried out piecewise in smaller intervals, allowing us to express the total evolution operator as the product of small time operators:

(10) |

where

(11) |

Evolution of the time dependent system is then reduced to the problem of approximating the propagator . Two approximations exist in the definition of , firstly that of approximating the matrix exponential and secondly the exact form of the matrix for which we wish to calulate the exponential. There are several methods for calulating the exponential of a matrix Moore (2011), here we use the simplest approximation, a Taylor expansion:

(12) |

Similarly there are many different approaches for deciding which matrix exponential to use as a propagator. Three approximations have been implemented: the so called exponential-midpoint propagator (EM), the enforced time-reversal symmetry (ETRS) propagator and the fourth order Magnus (M4) propagators, all of which are taken from the work of Marques et al.Castro, Marques, and Rubio (2004) on RT-TDDFT propagators, and are briefly described in our non-orthogonal basis for completeness.

The exponential midpoint propagator approximates the U by the exponential taken at :

(13) |

Implicitly enforcing time-reversibility, such that propagating forward from and backwards from by produce the same result, provides the so called enforced time-reversal symmetry method:

(14) |

Using the Magnus operator the exponential solution to Schrödinger equation for a time-dependent Hamiltonian may be written asMagnus (1954):

(15) |

where M is an infinite series of integrals providing an exact solution. Truncating this expansion to fourth order and approximating the integrals using Gauss-Legendre points as in Castro, Marques, and Rubio (2004) gives in our non-orthogonal basis:

(16) |

where .

It is important to note the presence of the inverse overlap matrix in these propagators, and again consider that while the overlap matrix will be well-ranged and suitable for truncation, the inverse overlap is not necessarily so. We therefore need to carefully test the sparsity of the product , and its effect on the unitarity of our propagators.

### ii.2 Linear Response

The idea behind extracting optical transitions from the linear response of a system to an external electric field is well knownYabana and Bertsch (1996); Tsolakidis, Sánchez-Portal, and Martin (2002). Propagating in real time provides direct access to the time-dependent charge density, and therefore the electronic response to external fields. Applying a time dependent external electric field polarised along axis j,

allows us to examine the time-dependent response of the system. Application of this electric field will produce an induced time-dependent dipole moment:

(17) |

As an example of the calculated repsonse of a system to an applied electric field, figure 1 illustrates the induced dipole response of a benzene molecule on application of a field with a Gaussian time profile, centred at .

Access to the time-dependent dipole moment allows us to calculate the time dependent polarisability:

The imaginary part of the polarisability is directly proportional to the absorption cross section, and the experimentally observed strength function, .

(18) |

As noted by Tsolakidis et. al., the approach satisifies the f-sum rule and the integration of the strength function over energy gives the number of electrons, which may be used as a measure of the completeness of the basis setTsolakidis, Sánchez-Portal, and Martin (2002).

Density matrix RT-TDDFT therefore has the potential to be an extremely useful tool for theoretically 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 several systems for which the electronic transitions have been studied experimentally and theoretically elsewhere, allowing us to make direct comparisons. For this purpose we have chosen four small molecules (Carbon monoxide, Methane, Ethylene and Benzene) and used our implementation to calculate the optical absorption spectra within the TDLDA approximation.

Meaningful comparison of our results with experiment requires the identification the electronic transitions to which the peaks in our calculated absorption spectra correspond. As we have mentioned in Casida’s approach information about electronic transitions is inherently produced, while in RT-TDDFT it is not.

It is often possible to identify the corresponding transition
by examining the polarisation and energy of peaks and comparing to
that of optically allowed transitions experimentally.
Where possible, in order to more confidently assign peaks of our calculated
absorption spectra to particular electronic transitions, we have followed the
procedure in Min, Cho, and Kim (2011) whereby a sinusoidal electric field tuned
to a particular excitation mode is applied. A resulting electronic resonance is set
up, allowing us to examine the difference between ground state charge
density and excited state charge density and thereby infer the electronic
transition.

### iii.1 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 Soler et al. (2002). These PAOs are eigenfunctions of the atomic pseudopotentials with a confinement energy shift used to determine a radial cut-off for the orbitals, beyond which they are zero. This confinement energy provides a single parameter to define the cut off radii for different orbitals, and is the energy 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 reproduced will generally not be satisfactory for calculating excited state properties, and therefore we illustrate the basis set dependence of two selected transitions for the CH molecule.

Multiple orbitals per angular momentum channel can be used (multiple-), with the shape of multiple orbitals determined by a split norm procedureSoler et al. (2002). This procedure uses a parameter to define the norm of a numerical orbital outside some radius where they match the tail of the first zeta PAO, and within this radius the vary smoothly to the origin. Subtracting this numerical orbital from the original PAO gives the multiple-zeta orbital. Of course it is possible to define these radii by hand and fine tune the basis set. In addition to multiple zeta, polarisation orbitals can be included within the basis set, and are obtained by solving the 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 so on. 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 calculations on the ethylene molecule with varying numbers of PAOs and two different confinement energies. The basis sets have been generated with a confinement energy of 1 meV and 5 meV, resulting in confinement radii of 4.93 and 4.24 Å for the carbon atoms respectively, and radii of 4.77 and 4.21 Å for the hydrogen atoms respectively. The total run time was 14.51 fs. (600 au.) with a time step of 0.0242 fs (0.1 au). The results can be seen in table 1.

Calculated energies for the transition show a wide variation with basis set choice, while the valence transition varies less. This is in line with expectation, given the more diffuse nature of the Rydberg transition we would expect its description to require larger basis. The effect of systematically increasing the number of basis functions is to improve our results with respect to that of the reference values. Similarly increasing the cut-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 increasing the range, will maximise the variational degrees of freedom available to describe our time dependent density matrix.

However our values are still far from those computed elsewhere, and we find generally that for small molecules it is essential to use a large basis with multiple extended polarisation orbitals in order to produce results in line with other works. In addition we find that fine tuning the radial cut-offs by hand, as opposed to using the confinement energy and split norm procedure, can allow us to improve the quality of our results for small molecules.

### iii.2 Small Molecule Results

Exhibited in table 2 are the calculated transitions for our four test molecules. In the case of the smallest molecules (carbon monoxide, ethylene, and methane) a hand tuned 5Z4P basis set is employed, while for benzene the result is obtained using a 2Z2P basis with a 5meV confinement energy (all the calculations satisfy the f-sum rule to ). Also presented in figure 2 are the optical absorption spectra for the benzene and carbon-monoxide molecules, along with the experimental data.

We can see a strong agreement between our results and that of other studies, giving us confidence in our implementation. Very good agreement is exhibited between the calculated benzene absorption spectra and the experimental values using a reasonably modest 2Z2P basis set. This highlights the point that for larger molecules we have generally found that the need for large hand tuned basis sets, as is necessary for the smaller molecules, is reduced. Typically results in agreement with those in the literature and experiment are found using smaller basis sets, a point that is important to bear in mind, given the context of linear scaling methods.

## Iv Propagator Unitarity

Having demonstrated the correctness of our implementation and explored the influence of basis sets, we now turn to our main concern, the effects of localisation in linear scaling methods on the accuracy of results.

We wish the total charge in our system to remain stable, and in order for this to be the case the propagators must be unitary with respect to the non-orthogonal basis set:

(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 it were exact, our propagators would indeed exhibit this property. However, as it is impossible for us to store an infinite sum on our computer, we must truncate our Taylor expansion at some point. Doing so will introduce errors, with two factors affecting the scale of the break from unitarity; the time step and the number of terms in our summation. While we can extend our expansion arbitrarily, and reduce the time step arbitrarily, we wish to avoid excess computational expense by keeping the expansion as small as possible and the time step as large as possible within some acceptable margin of accuracy. We can directly examine the unitarity of our propagators through equation 19.

### iv.1 Time-Step Dependence

As a test we have examined the extent of the break from unitarity for a range of time-steps and number of terms in the matrix exponential expansion. We have used a small molecule for the purpose, benzene, with a small applied electric field perturbation with a Gaussian profile centered on .

Exhibited in figure 3 we can see the dependence on simulation time step of the propagator unitarity, with the obvious trend being that as the time step is reduced the propagator approaches unitarity. We can see that even for time steps up to a.u. the propagator maintains its unitarity to a high degree (similar results were obtained for each of the propagators). 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 deteriorates with the propagation eventually becoming unstable for large timesteps. While the maximum permissible timestep will depend on the system under study, we found that generally a timestep of 0.06 a.u. or below provided satisfactory charge conservation.

The form of our propagators requires the extrapolation of the Hamiltonian matrix to some unknown point beyond the current time t, . As suggested by Marques et al.Castro, Marques, and Rubio (2004) in order to minimise errors it is possible to carry this procedure out self-consistently. In our case meaning that we propagate K to K based on an extrapolated Hamiltonian. We then construct a new Hamiltonian matrix H using K. can then be interpolated from Hamiltonian matrices for times up to and including , and the whole procedure is iterated until some self-consistency criteria is obtained. Generally speaking this procedure 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-consistency procedure on the charge conservation can be seen in figure 4. While the self-consistency cycle is found to improve the charge conservation, in reality for small time steps the difference in charge conservation and calculated properties is not found to be significant enough to warrant the extra computational load of constructing the Hamiltonian matrix several times per time-step. As a compromise we enforce the self-consistency only for a small number of steps () at the beginning of a run, typically when our external electric field is applied for the study of the linear response and the external perturbation is largest.

A significant point to note is that little difference is exhibited between the calculated results using each of the three propagators in terms of charge conservation, and in general we have found this to be the case. It is reported that for systems with strongly time-dependent Hamiltonians the fourth order Magnus propagator, U, is advantageousCastro, Marques, and Rubio (2004), but for our present work this is not the case and we have opted for the simplest exponential midpoint propagator throughout.

### iv.2 Matrix Exponential Truncation

The effect of truncating the Taylor expansion used to evaluate the matrix exponential on the unitarity of the propagator can be seen in figure 6. We see that reducing the number of terms reduces the unitarity of the propagator, as expected. Looking at figure 5 the convergence of the charge conservation with the number of terms in the exponential expansion can 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.

## V Alkane Molecules: Testing Matrix Truncation Effects

In this section we perform calculations on long chain alkane molecules.Our aim is to examine the effect of matrix truncation on the propagation of the density matrix and propagator unitarity, along with the computational scaling with system size.

As a first step we calculate the absorption spectra for the CH molecule for several different basis sets using the generalised gradient PBE functionalPerdew, Burke, and Ernzerhof (1996) (all further calculations in 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 is found to reduce, and the reported adsorption onset for CH is 175 nm. Costner et al. (2009) ( 7.1 eV). We see that as the number of PAOs in the basis set is increased the calculated absorption onset approaches this value. Particularly noticeable is the change of the absorption energy caused by the addition of polarisation orbitals. Similarly a significant shift is induced by extending the range of the PAOs (a variation from 55 meV to 25 meV in the confinement energy extends the radii of the carbon and hydrogen basis sets by 0.35 Å and 0.33 Å respectively). This is understandable, given that the first transitions in the alkane molecules are reported as being Rydberg in characterCostner et al. (2009), we would expect the addition of more diffuse PAOs to improve the description of these excitations. Given the well documented difficulties of TDDFT to accurately describe Rydberg transitions Casida et al. (1998), and given that this is not our aim in any case, we proceed to carry out our tests with the SZP and SZ2P basis sets generated using a confinement energy of 55meV (radial cut off for the PAOs is 3.31Å and 3.12Å for carbon and hydrogen respectively).

(i) | (ii) |

Yam et al. have previously studied the long chain alkanes within the linear scaling excited state regime Yam, Yokojima, and Chen (2003b), calculating the absorption onset at around 8 eV for CH with the LDA functional. However little discussion of the effects of matrix truncation on propagator unitarity have been presented elsewhere.

### v.1 Propagator Truncation

The use of a basis of non-orthogonal atomic orbitals requires the inverse overlap matrix for our propagation (indeed this matrix is required for ground state calculations in any case), as seen in equation 11. In order to compute the inverse overlap matrix Conquest uses Hotelling’s methodPalser and Manolopoulos (1998), however for poorly conditioned overlap matrices computing the inverse overlap matrix can prove difficult. In our current implementation of TDDFT the atoms remain stationary and so too, therefore, does the overlap matrix. Therefore we have also included the possibility of computing the inverse overlap with the SCALAPACK routines. Although the scaling will not be linear, computing the inverse overlap in this way makes only a relatively small contribution to our total TDDFT runtime, as we only calculate the inverse overlap once at .

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 tested the effect of truncating both the S matrix and the SH matrix on the propagation. Figure 8 shows the average absolute error in the matrix elements of S and the SH matrices caused by truncation (the error in S elements given is the average of the elements of the SS-I matrix, and the error in the SH is calculated with the values from an untruncated S matrix).

As the range of the matrices increases the error caused by the truncation converges towards zero, as we expect. The S matrix converges less quickly than the SH matrix, indicating that it is more dense than the SH matrix. The effect truncation of these matrices has on the unitarity of the propagators can be seen in figure 9. We see that the unitarity converges as the SH range increases, and the propagators are converged with a range of around 22.5-27.5 Bohr. This indicates that the SH matrix is indeed sparse, while the S matrix is less so, and we can safely truncate it. It is important to not that we don’t explicitly use the S in our propagators, only the SH matrix. Although it makes sense to truncate the S matrix, given that we are truncating SH and that the Hamiltonian matrix is sparse. We can see this by noting that the unitarity of the propagator in figure 9 is also well converged for each of the truncation ranges imposed on the inverse overlap.

As additional atoms are added the Hamiltonian matrix, overlap matrix, and the inverse overlap will vary. Increasing the system size may therefore affect the ranges of these matrices. While we only use the SH matrix in our calculation, comparison of the density of both matrices have been included. We have tested this effect by fixing the S and SH ranges at 30 and 35 Bohr respectively, and examined the error in the truncated SH matrix with system size with the results shown in figure 10. We see that the error changes slightly on increasing system size, but converges as the size increases. Consequently the propagator unitarity was found to exhibit the same trend. This illustrates that the SH 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 made to 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 of this effect we have examined the CH molecule with a larger basis set (SZ2P as opposed to SZP). Exhibited in figure 11 is the absolute value of the matrix with SH matrix truncation range. Despite the larger number of basis set functions we see that the SH matrix is still well ranged, although the range is wider when compared to the SZP results of figure 9, and again a truncation will lead to a computational load that scales linearly with system size.

A further point to note is that it is possible to avoid the use of the inverse overlap matrix in the TDDFT propagation altogether. Yam et al. have employed a Cholesky orthogonalisation scheme to bypass the need for the inverse overlapYam, Yokojima, and Chen (2003b).n However using this scheme requires the inverse of the Cholesky decomposition, and it is not apparent that it will be more sparse than the inverse overlap. It is possible that this scheme might improve the calculation of the propagator, as the orthogonalised Hamiltonian may be more localised than our SH matrix. Calculating the Cholesky decomposition can be made to scale linearly, and implementation of this alternative method has already begun in order to contrast the two approaches. However the parallelisation of Cholesky inversion is difficult given the Conquest matrix storage, and inversion of the overlap matrix remains important.

### v.2 Density Matrix Truncation, Scaling and Limits

Finally we examine the effect of truncating the density matrix, and have performed calculations generating spectra for the CH molecule at varying truncation radii, R, of the density matrix. Typically for ground state calculations a suitable typical density matrix truncation range is around 16-20 Bohr. The results can be seen in figure 12, and generally 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 from the comparison of and that there is good agreement for the initial transitions, as well as the general shape of the spectra.

Applying this Bohr cut-off (along with a cut off of 35 Bohr. on the matrix) we can examine the computational scaling with system size, with the results exhibited in figure 13. Clear linear scaling of the computational workload up to well over 1000 atoms is exhibited, illustrating the potential power of the method.

Finally a few comments on the limits of the approach must be made. TDDFT for long-range charge transfer is well known to be poorly described by local and semi-local functionalsDreuw, Weisman, and Head-Gordon (2003). While we have employed LDA and GGA functionals here, linear scaling exact exchange has also been recently implemented in the Conquest code, allowing the use of non-local functionals with this approach in the future.

While the near-sightedness principle dictates that the ground-state density matrix is exponentially localised for well gapped systems, there is no formal justification for the localisation of the response density matrix. As noted in Zuehlsdorff et al. (2013), for systems with well localised excitations it would be expected that the response density matrix could be truncated safely and linear scaling achieved, while for systems with delocalised excitations this will not be the case.

## Vi Conclusions

We have outlined our implementation of real-time time dependent density functional theory in the Conquest (N) code. We have demonstrated the soundness of the implementation through benchmark tests for small molecules, and also discussed the effect of basis set and system sizes on the results.

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 size by its truncation. We have discussed the range of our propagator matrices for an alkane chain test system, and the implications of this matrix truncation on the unitarity of the propagation. Similarly we have examined the effect of truncating the density matrix on the calculated optical absorption spectra, showing that the range required is much more extended than that required for converged ground state properties. Nevertheless, we have shown that accurate linear scaling TDDFT calculations are practical. While the impact of localisation cut-off in the charge density matrix on these TDDFT calculations is a topic warranting further study, we have shown that in truncating these matrices at a suitable point we obtain a computational load that increases linearly with system size. This offers a complementary approach to the usual Casida linear response approach: linear response TDDFT is well suited to relatively small systems, while linear scaling RT-TDDFT offers a viable method for studying excitations in large systems. We have shown linear scaling beyond 1,000 atoms, and 10,000+ atoms are perfectly practical with the excellent parallel scaling available in Conquest.

## Acknowledgements

C.O’R. is supported by the MANA-WPI project and D.R.B. was funded by the Royal Society. We thank Umberto Terranova for useful discussions. This work made use of the facilities of ARCHER, the UK’s national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and Technology through EPSRC’s High End Computing Programme. Calculations were performed at HECToR through the UKCP Consortium. The authors acknowledge the use of the UCL Legion High Performance Computing Facility, and associated support services, in the completion of this work.

## References

- Bowler and Miyazaki (2012) D. R. Bowler and T. Miyazaki, “(n) methods in electronic structure calculations,” Reports on Progress in Physics 75, 036503 (2012).
- Bowler and Miyazaki (2010) D. R. Bowler and T. Miyazaki, “Calculations for millions of atoms with density functional theory: linear scaling shows its potential,” Journal of Physics: Condensed Matter 22, 074207 (2010).
- VandeVondele, Bor¿tnik, and Hutter (2012) J. VandeVondele, U. Bor¿tnik, and J. Hutter, “Linear scaling self-consistent field calculations with millions of atoms in the condensed phase,” Journal of Chemical Theory and Computation 8, 3565–3573 (2012), http://dx.doi.org/10.1021/ct200897x .
- Casida (1996) M. E. Casida, in Recent Developments and Applications of Modern Density Functional Theory (Elsevier, Amsterdam, 1996).
- Kronik et al. (2006) 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 electronic structure calculations: recent advances and novel applications to nano-structures,” physica status solidi (b) 243, 1063–1079 (2006).
- Rocca et al. (2008) D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, “Turbo charging time-dependent density-functional theory with lanczos chains,” The Journal of Chemical Physics 128, 154105 (2008).
- Zuehlsdorff et al. (2013) 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 response formalism,” The Journal of Chemical Physics 139, 064104 (2013).
- Yabana and Bertsch (1996) K. Yabana and G. F. Bertsch, ‘‘Time-dependent local-density approximation in real time,” Phys. Rev. B 54, 4484–4487 (1996).
- Yokojima and Chen (1998) S. Yokojima and G. Chen, “Time domain localized-density-matrix method,” Chemical Physics Letters 292, 379 – 383 (1998).
- Yokojima and Chen (1999a) S. Yokojima and G. Chen, “Linear scaling calculation of excited-state properties of polyacetylene,” Phys. Rev. B 59, 7259–7262 (1999a).
- Yokojima and Chen (1999b) S. Yokojima and G. Chen, “Linear-scaling localized-density-matrix method for the ground and excited states of one-dimensional molecular systems,” Chemical Physics Letters 300, 540 – 544 (1999b).
- Yokojima, Zhou, and Chen (1999) S. Yokojima, D. Zhou, and G. Chen, “Linear-scaling computation of ground state with time-domain localized-density-matrix method,” Chemical Physics Letters 302, 495 – 498 (1999).
- Liang, Yokojima, and Chen (1999) W. Liang, S. Yokojima, and G. Chen, “Generalized linear-scaling localized-density-matrix method,” The Journal of Chemical Physics 110, 1844–1855 (1999).
- Liang et al. (2000) W. Liang, S. Yokojima, D. Zhou, and G. Chen, “Localized-density-matrix method and its application to carbon nanotubes,” The Journal of Physical Chemistry A 104, 2445–2453 (2000), http://pubs.acs.org/doi/pdf/10.1021/jp990818a .
- Tsolakidis, Sánchez-Portal, and Martin (2002) A. Tsolakidis, D. Sánchez-Portal, and R. M. Martin, “Calculation of the optical response of atomic clusters using time-dependent density functional theory and local orbitals,” Phys. Rev. B 66, 235416 (2002).
- Yam et al. (2012) C. Yam, Q. Zhang, F. Wang, and G. Chen, “Linear-scaling quantum mechanical methods for excited states,” Chem. Soc. Rev. 41, 3821–3838 (2012).
- Yam, Yokojima, and Chen (2003a) C. Y. Yam, S. Yokojima, and G. Chen, “Localized-density-matrix implementation of time-dependent density-functional theory,” The Journal of Chemical Physics 119, 8794–8803 (2003a).
- Yam, Yokojima, and Chen (2003b) C. Yam, S. Yokojima, and G. Chen, “Linear-scaling time-dependent density-functional theory,” Phys. Rev. B 68, 153105 (2003b).
- Ratcliff, Hine, and Haynes (2011) L. E. Ratcliff, N. D. M. Hine, and P. D. Haynes, “Calculating optical absorption spectra for large systems using linear-scaling density functional theory,” Phys. Rev. B 84, 165131 (2011).
- Goringe et al. (1997) C. Goringe, E. Hernández, M. Gillan, and I. Bush, “Linear-scaling dft-pseudopotential calculations on parallel computers,” Computer Physics Communications 102, 1 – 16 (1997).
- Moore (2011) G. Moore, “Orthogonal polynomial expansions for the matrix exponential,” Linear Algebra and its Applications 435, 537 – 559 (2011), ¡ce:title¿Special Issue: Dedication to Pete Stewart on the occasion of his 70th birthday¡/ce:title¿.
- Castro, Marques, and Rubio (2004) A. Castro, M. A. L. Marques, and A. Rubio, “Propagators for the time-dependent kohn–sham equations,” The Journal of Chemical Physics 121, 3425–3433 (2004).
- Magnus (1954) W. Magnus, “On the exponential solution of differential equations for a linear operator,” Communications on Pure and Applied Mathematics 7, 649–673 (1954).
- Matsuzawa et al. (2001a) N. N. Matsuzawa, A. Ishitani, D. A. Dixon, and T. Uda, “Time-dependent density functional theory calculations of photoabsorption spectra in the vacuum ultraviolet region,” The Journal of Physical Chemistry A 105, 4953–4962 (2001a), http://pubs.acs.org/doi/pdf/10.1021/jp003937v .
- Hammond and Price (1955) V. J. Hammond and W. C. Price, “Oscillator strengths of the vacuum ultra-violet absorption bands of benzene and ethylene,” Trans. Faraday Soc. 51, 605–610 (1955).
- Min, Cho, and Kim (2011) S. K. Min, Y. Cho, and K. S. Kim, “Efficient electron dynamics with the planewave-based real-time time-dependent density functional theory: Absorption spectra, vibronic electronic spectra, and coupled electron-nucleus dynamics,” The Journal of Chemical Physics 135, 244112 (2011).
- Hu, Sugino, and Miyamoto (2006) 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 74, 032508 (2006).
- Chan, Cooper, and Brion (1993) W. Chan, G. Cooper, and C. Brion, “Absolute optical oscillator strengths for discrete and continuum photoabsorption of carbon monoxide (7-200 ev) and transition moments for the system,” Chemical Physics 170, 123 – 138 (1993).
- Matsuzawa et al. (2001b) N. N. Matsuzawa, A. Ishitani, D. A. Dixon, and T. Uda, ‘‘Time-dependent density functional theory calculations of photoabsorption spectra in the vacuum ultraviolet region,” The Journal of Physical Chemistry A 105, 4953–4962 (2001b), http://pubs.acs.org/doi/pdf/10.1021/jp003937v .
- Curtis and Walker (1989) M. G. Curtis and I. C. Walker, “Low-energy electron-impact excitation of methane, silane, tetrafluoromethane and tetrafluorosilane,” J. Chem. Soc., Faraday Trans. 2 85, 659–670 (1989).
- Yabana and Bertsch (1999) K. Yabana and G. F. Bertsch, “Time-dependent local-density approximation in real time: Application to conjugated molecules,” International Journal of Quantum Chemistry 75, 55–66 (1999).
- Koch and Otto (1972) E. Koch and A. Otto, “Optical absorption of benzene vapour for photon energies from 6 ev to 35 ev,” Chemical Physics Letters 12, 476 – 480 (1972).
- Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, ‘‘The siesta method for ab initio order- n materials simulation,” Journal of Physics: Condensed Matter 14, 2745 (2002).
- Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
- Costner et al. (2009) 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 cyclic alkanes: Vuv absorbance and index of refraction,” The Journal of Physical Chemistry A 113, 9337–9347 (2009), pMID: 19630422, http://pubs.acs.org/doi/pdf/10.1021/jp903435c .
- Casida et al. (1998) M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, ‘‘Molecular excitation energies to high-lying bound states from time-dependent density-functional response theory: Characterization and correction of the time-dependent local density approximation ionization threshold,” The Journal of Chemical Physics 108, 4439–4449 (1998).
- Palser and Manolopoulos (1998) A. H. R. Palser and D. E. Manolopoulos, “Canonical purification of the density matrix in electronic-structure theory,” Phys. Rev. B 58, 12704–12711 (1998).
- Dreuw, Weisman, and Head-Gordon (2003) A. Dreuw, J. L. Weisman, and M. Head-Gordon, “Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange,” The Journal of Chemical Physics 119, 2943–2946 (2003).