All-electron full-potential implementation of real-time TDDFT in exciting
aa r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Submitted to:
Electronic Structure
All-electron full-potential implementation ofreal-time TDDFT in exciting
Ronaldo Rodrigues Pela , , Claudia Draxl , Physics Department and IRIS Adlershof, Humboldt-Universit¨at zu Berlin, ZumGroßen Windkanal 2, 12489 Berlin European Theoretical Spectroscopy Facility (ETSF)E-mail: [email protected]
Abstract.
Linearized augmented planewaves combined with local-orbitals (LAPW+lo)are arguably the most precise basis set to represent Kohn-Sham states. When em-ployed within real-time time-dependent density functional theory (RT-TDDFT), theypromise ultimate precision achievable for exploring the evolution of electronic excita-tions. In this work, we present an implementation of RT-TDDFT in the full-potentialLAPW+lo code exciting . We benchmark our results against those obtained by linear-response TDDFT with exciting and by RT-TDDFT calculations with the Octopuscode, finding a satisfactory level of agreement. To illustrate possible applications ofour implementation, we have chosen three examples: the dynamic behavior of excita-tions in MoS induced by a laser pulse, the third harmonic generation in silicon, anda pump-probe experiment in diamond. ll-electron full-potential implementation of real-time TDDFT in exciting
1. Introduction
TDDFT is a powerful tool to study excitations in many-body systems [1]. It is formallyan exact theory that guarantees the existence of a one-to-one mapping between theevolution of the many-body wavefunction and a much less complicated object, namelythe time-dependent density [1, 2]. Compared to many-body perturbation theory basedon Green functions, TDDFT is computationally less demanding, which allows forstudying systems containing up to hundreds or even thousands of atoms [3, 4, 5, 6, 7].In practical problems, TDDFT is usually employed either in the linear-response(LR) regime, where the density is evaluated in the frequency domain as a first-orderresponse to an external perturbation potential, or directly in the time domain, byevolving the Kohn-Sham (KS) wavefunctions [5]. Both approaches have advantagesand limitations. Here, we focus on RT-TDDFT that allows, among others, forassessing the nonlinear regime, for studying the dynamics of excitations in responseto ultra-fast laser pulses as observed in pump-probe spectroscopy, and for observingthe coupling of electronic excitations to the vibrations of the nuclei in real time[8, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16].On the one hand, a key point for the accuracy of TDDFT and, thus a mostmeaningful comparison with experiments, is the choice of the time-dependent exchange-correlation (TD-XC) functional. Therefore, the improvement of existing approximationshas been in the focus of recent active research [5, 17, 18, 19, 20, 21]. On the otherhand, for a given XC functional, it is desirable to achieve ultimate numerical precisionwithin a calculation. This usually translates into a question about the quality of thebasis employed to represent KS wavefunctions. In this sense, with its all-electron, full-potential LAPW+lo approach, exciting has proven to be one of the most accurate ab initio codes [22], capable of reaching micro-Hartree precision [23]. Furthermore, itis a user-friendly code with a growing community, and offers plenty of tutorials on itsimplementations [24]. Nevertheless, up to now, TDDFT is available in exciting onlyin its LR formulation [25] while RT-TDDFT is still missing.In this paper, we fill this gap and present our implementation of RT-TDDFT in the exciting code. A summary of the theory behind and the details of the implementationare given in Section 2. In Section 3, we provide benchmarks by comparing our resultswith those obtained by the Octopus code [26], while more examples can be foundin the Appendices. In Section 4, we demonstrate three interesting features of ourimplementation, namely (1) an analysis of the behavior of MoS after excitation with alaser pulse, (2) a study of the third harmonic generation in silicon, and (3) a simulation ofa pump-probe experiment in diamond. Finally, in Section 5, we provide our conclusions.
2. Theory and Implementation
We start by considering a physical system with periodic boundary conditions subjectedto an electric field E ( t ) with spatial variations on a scale much larger than the periodicity, ll-electron full-potential implementation of real-time TDDFT in exciting E ( t ) in the KSHamiltonian would be by addition of a dipole term r · E ( t ) that would break the desiredperiodicity. In this case, it is advantageous to employ the velocity-gauge [27, 28, 9, 29]:ˆ H ( r , t ) = 12 (cid:18) − i ∇ + 1 c A ( t ) (cid:19) + v KS ( r , t ) , (1)where A ( t ) is the vector potential, given in this gauge by A ( t ) = − c R t E ( t ′ ) dt ′ , c is thespeed of light, and v KS is the TD-KS potential, a sum of the TD ionic, Hartree andXC potentials. We assume here the adiabatic approximation for the TD-XC potential[30, 31]. A KS wavefunction ψ j k ( r , t ) labeled with index j and associated to a wavevector k evolves as ˆ H ( r , t ) ψ j k ( r , t ) = i ∂∂t ψ j k ( r , t ) . (2)In exciting , each KS wavefunction is expanded in terms of the LAPW+lo basisset with coefficients C j k ( t ): | ψ j k ( t ) i = X G C j kG ( t ) | φ G + k i + X γ C j k γ ( t ) | φ γ i , (3)where | φ G + k i and | φ γ i represent the LAPW part of the basis and the local orbitals(lo), respectively; G is a reciprocal lattice vector. With the definition of the basis,the integration of Eq. (2) reduces to the problem of finding how C j k ( t ) propagate intime. For this, many approaches are available [32, 33, 34]. Apart from the classicalRunge-Kutta method for differential equations, we have implemented the followingpropagators: (i) simple exponential, (ii) exponential at the midpoint, (iii) approximateenforced time-reversal symmetry, (iv) commutator-free Magnus expansion of 4th order,and (v) exponential using a basis of the Hamiltonian’s eigenvectors. To illustrate how ourimplementation works, we choose the simple exponential propagator, while Appendix Bdetails the other cases.The evolution of KS wavefunctions in terms of a propagator ˆ U ( t + ∆ t, t ) is: | ψ j k ( t + ∆ t ) i = ˆ U ( t + ∆ t, t ) | ψ j k ( t ) i , (4)where ˆ U ( t + ∆ t, t ) = ˆ T (cid:20) exp (cid:18) − i Z t +∆ tt d τ ˆ H ( τ ) (cid:19)(cid:21) , (5)ˆ T being the time-ordering operator. For the simple exponential propagator, ˆ U ( t + ∆ t, t )is regarded purely as exp[ − i∆ t ˆ H ( t )], the approximation being better, the smaller thetime step ∆ t is. For this propagator, the following expression dictates the evolution of C j k ( t ): C j k ( t + ∆ t ) = exp[ − i∆ tS − k H k ( t )] C j k ( t ) . (6)The matrix exponential in Eq. (6) is approximated by a Taylor expansion up to theorder defined by the user (4 is the default). H k and S k are, respectively, the Hamiltonianand overlap matrices in the basis set. Since our basis consists of two parts, i.e. , LAPWs ll-electron full-potential implementation of real-time TDDFT in exciting H k and S k usually have a block structure, exemplified for the Hamiltonian, asfollows: " h φ k + G | ˆ H ( t ) | φ k + G ′ i h φ k + G | ˆ H ( t ) | φ γ ′ ih φ γ | ˆ H ( t ) | φ k + G ′ i h φ γ | ˆ H ( t ) | φ γ ′ i . (7)In exciting , each block is calculated employing different strategies, as described in Ref.[25]. If µ and ν denote generic indexes that can be associated to a LAPW or lo, then,following Eq. (1), an arbitrary element [ H k ( t )] µν = h φ µ | H ( t ) | φ ν i can be written as:[ H k ( t )] µν = 12 h∇ φ µ |∇ φ ν i + h φ µ | v KS ( t ) | φ ν i ++ A ( t )2 c [ S k ] µν − i c A ( t ) · h φ µ |∇| φ ν i . (8)The procedure to obtain h∇ φ µ |∇ φ ν i , h φ µ | v KS ( t ) | φ ν i , and [ S k ] µν is detailed in Ref. [25],whereas the momentum matrix elements, h φ µ | − i ∇| φ ν i , are calculated as described inRef. [35]. The dielectric and optical properties are quantities that can be measured by variousexperimental probes. They are also of main interest in LR-TDDFT calculations. Toobtain them with RT-TDDFT, the behavior of the macroscopic current density, J ( t ), inresponse to an external field needs to be evaluated. For the case of local and semilocalKS functionals, J ( t ) can be obtained as J ( t ) = iΩ X j k w k f j k h ψ j k ( t ) |∇| ψ j k ( t ) i − N A ( t ) c Ω , (9)where N is the number of valence electrons in the unit cell with volume Ω, w k is theweight of the considered k -point, and f j k is the occupation number of the correspondingKS state. After Fourier transform, we can obtain the components of the opticalconductivity σ and the dielectric tensor ε as σ αβ ( ω ) = J α ( ω ) E β ( ω ) , ε αβ ( ω ) = δ αβ + 4 π i σ αβ ( ω ) ω , (10)where the indexes α , β mean the cartesian directions x , y or z . It is convenient toconsider an impulsive electric field E δ ( t ) in a specific direction α such that E α ( ω ) = E .This is known as the transverse geometry [9]. To understand its physical meaning, welook at the interface between the studied material and the vacuum, where the electricfield comes from. Transverse geometry means that the field direction is parallel to theinterface. Conversely, in the longitudinal geometry the electric field is perpendicularto the interface [9]. In this case, the displacement fields D ( t ) inside and outside thesystem are related to each other through the surface charge, as given by the boundaryconditions for electromagnetic fields [36]. Following Refs. [9], [27], and [28], we considerthe external component of the vector potential in the longitudinal geometry as given ll-electron full-potential implementation of real-time TDDFT in exciting A ext ( t ) = − cD θ ( t ) e α or, equivalently, by D ( t ) = − (1 /c ) d A ext /dt = D δ ( t ) e α . Theinduced vector potential, A ind ( t ), is obtained from the current density as d A ind dt = 4 πc J ( t ) , (11)and the total vector potential A ( t ), appearing in the Hamiltonian, is calculated as thesum A ext ( t ) + A ind ( t ). A quantity of interest is the number of excited electrons after the interaction with alaser pulse [37, 11, 38, 39, 9, 40, 41, 42]. In RT-TDDFT, the occupation number f j k of a KS state is kept fixed to its initial value. As the wavefunctions evolve, they areno longer eigenstates of ˆ H ( t ). It is possible to describe the number of excitations byprojecting | ψ i k ( t ) i onto the adiabatic ground state of ˆ H ( t )’s eigenfunctions [42, 11, 38]or onto the reference ground state at t = 0 [42, 9, 41]. We opt here for the latter.Therefore, for a given k -point, we define the number of electrons that have been excitedto an unoccupied KS state, labeled j , as m ej k ( t ) = X i f i k |h ψ j k (0) | ψ i k ( t ) i| . (12)Similarly, the number of holes created in an occupied KS state j ′ can be specified as m hj ′ k ( t ) = f j ′ k − X i f i k |h ψ j ′ k (0) | ψ i k ( t ) i| . (13)Thus, the total number of excited electrons in a unit cell can be obtained by consideringall the unoccupied states N exc ( t ) = j unocc X j k w k m ej k ( t ) = j ′ occ X j ′ k w k m hj ′ k ( t ) . (14) We follow the same parallelization strategy as already adopted in other parts of exciting , i.e. , over k -points [25]. In Fig. 1, we contrast the performance of twodifferent levels of parallelization for calculations carried out on a single node withmultiple processors: Open Multi-Processing (OpenMP) and Message Passing Interface(MPI). Although the speedup in both cases appears to be very close to the ideal one,MPI alone tends to be more efficient – also when compared to a hybrid parallelization(using both OpenMP and MPI, not shown in the figure). In the inset of Fig. 1, wedepict the speedup of MPI when the calculations are distributed over a higher numberof nodes, still showing fairly close to ideal scaling (speedup of 187 for 256 processors). ll-electron full-potential implementation of real-time TDDFT in exciting 𝘕 𝘵𝘩𝘳𝘦𝘢𝘥𝘴 or 𝘕 𝘱𝘳𝘰𝘤 Sp ee d u p OpenMPMPIIdeal100 200100200
Figure 1.
Comparison between the two parallelization schemes MPI and OpenMPfor the current density in diamond using 16 × ×
16 irreducible k -points. The insetdepicts the speedup by MPI when the job is distributed among several nodes. Forcomparison, the ideal scaling is indicated by the dashed line. In the following, we analyze the impact of the three most important parametersgoverning the precision of RT-TDDFT calculations, namely: the time step, the numberof k -points, and the size of the basis. For a given value p of any of these parameters,we adopt the root-mean square error (RMSE) E p = s R T ( j p ( t ) − j ref ( t )) d tT , (15)to address the convergence behavior. Here, j ref is a reference value for the currentdensity, corresponding to the optimal parameter, and T stands for the end-time, up towhich the evolution of KS wavefunctions is considered. We consider diamond under an impulsive displacement field along the[001], direction given by D = 0 . δ ( t −
1) in atomic units (a.u.). In Fig. 2, we show thecurrent-density response to this field in the same direction. On the one hand, the currentdensity is apparently insensitive to the time step, as the various curves seem to coincide,suggesting swift convergence of the results with decreasing time-step. Interestingly, timesteps above 0.2 a.u. lead to divergence. On the other hand, taking the current densityobtained with a time step of 0.001 a.u. as j ref in Eq. (15), the RMSE depicted on theright side of the figure shows that the calculations become indeed more precise whenthe time step is reduced. And there is actually no saturation behavior, i.e. , the RMSEscales with the time step by a power law of 0.8. Although the value of the exponentdepends on the material and on the method employed to propagate the wavefunctions,such power laws are found as a quite general trend, as already pointed out in Ref. [32].A similar conclusion can be drawn in the case of silicon exposed to an electric field of ll-electron full-potential implementation of real-time TDDFT in exciting Time [a.u.] −0.002−0.0010.0000.0010.002 C u rr e n t d e n s i t y [ a . u . ] Time step [a.u.]0.0010.005 0.010.02 0.050.1 10 −2 −1 Time step [a.u.] −4 −3 R M S E [ a . u . ] Figure 2.
Convergence behavior of the current density in diamond with respect tothe chosen time step. The right panel shows the RMSE, taking the calculation withthe smallest time step as reference. two different forms (see Figs. C1 and C2 in Appendix C).
Time [a.u.] −0.0010.0000.0010.002 C u rr e n d e n s i y [ a . u . ] Number of 𝗸-points20 Number of 𝘬-points −6 −5 −4 R M S E [ a . u . ] Figure 3.
Impact of the number of k -points on the convergence of the currentdensity. Results for diamond under the influence of an impulsive displacement field of D = 0 . δ ( t −
2) applied along the [001] direction. k -points Analogous to ground-state calculations, the k -grid has adirect impact on the quality of the current density (Eq. (9)) and the time-dependentelectronic density, calculated as n ( r , t ) = X j k w k f j k | ψ j k ( r , t ) | . (16)To illustrate its role, we consider diamond exposed to an impulsive displacement field D = 0 . δ ( t −
2) a.u. along the [001] direction. Figure 3 depicts how the numberof k -points affects the current density. Once more, the RMSE follows a power-lawdependence on the investigated parameter (now, the number of k -points), as seen in the ll-electron full-potential implementation of real-time TDDFT in exciting Figure 4.
Influence of the k -grid on the current density with (left) and without (right)an offset. Another relevant aspect concerning the k -grid is a possible offset that usually lowersthe symmetry, leading to a set of symmetrically inequivalent k -points. When the goal is, e.g. , to obtain the dielectric function, the offset helps to avoid symmetrically redundantcontributions. To exemplify the effect of such offset, we take as a test case siliconexposed to an electric field along the [001] direction given by E ( t ) = 0 . δ ( t − . k -grids with anoffset of 0 . b + 0 . b + 0 . b (where b i are the reciprocal lattice vectors). Thegraph on the right side shows the case when no offset is taken into account. Comparingboth graphs, we verify that, apart from a vertical shift, the current density convergesfaster with respect to the k -grid when an offset is considered. The vertical shift signalsthat the offset induces an artificial long-time behavior J ( t → ∞ ) that does not convergeto zero for coarser grids. This is not the case without an offset. Actually, in thelinear regime, the summation in (9) should be ideally zero when the excitation fieldis removed. An offset may erroneously hamper cancellation of terms, this effect beingmuch less pronounced when finer k -grids are considered.Figure 5 shows an equivalent comparison for the imaginary part of the dielectricfunction, calculated from the Fourier transform of the current density, as given in Eq.(10). We note a “fake” plasmon at smaller frequencies (0-2 eV). This has already beenreported in the literature as a consequence of the velocity-gauge [27, 29, 9, 43]. Whenno offset is included, the convergence with respect to the number of k -points is slower.In contrast, when an offset is taken into account, calculations with a k -grid of 8 × × k -points can describe those transitions thatoccur in the vicinity of k -points with high-symmetry. ll-electron full-potential implementation of real-time TDDFT in exciting Figure 5.
Imaginary part of the dielectric function of Si varying the number of k -points with (left) and without (right) an offset to break the symmetry. In the case of LAPW+lo, the dimensionless parameter rgkmax togetherwith the number of lo’s determine the quality of the basis. We need to inspect theirimpact on the convergence behavior separately. Starting with rgkmax , we consider thecurrent density in diamond exposed to an electric field E ( t ) = 0 . δ ( t −
1) a.u. along[001]. From Fig. 6, we conclude that the RMSE decreases exponentially when increasing rgkmax . A similar behavior can also be observed in Fig. C6 (Appendix C) for siliconexposed to a sinusoidal electric field modulated by a gaussian-like envelope.
Time [a.u.] −0.0010.0000.0010.002 C u rr e n t d e n s i t y [ a . u . ] rgkma 109.59.0 8.58.07.5 7.06.56.0 5.55.0 4.54.0 5.0 7.5 10.0 rgkma [a.u.] −6 −5 −4 −3 R M S E [ a . u . ] Figure 6.
Convergence behavior of the current density in diamond for different valuesof rgkmax , determining the basis-set size. The corresponding RMSE is displayed onthe right.
To check the role of lo’s, we consider silicon subjected to the electric field E ( t ) = E m sin (cid:20) π ( t − t ) T pulse (cid:21) cos( ω t ) , (17)along [001] for t ≤ t ≤ T pulse , and 0 otherwise. This function describes a periodicwave with angular frequency ω modulated by a gaussian-like function, correspondingto a laser shape frequently employed in experiment. We choose, E m = 1, ω = 0 . ll-electron full-potential implementation of real-time TDDFT in exciting t = 2, and T pulse = 452, all quantities in a.u. In Fig. 7, we show how the current densitychanges when enhancing the basis with more lo’s. Adding lo’s with p or d charactertends to improve the precision more than lo’s with s character. lo’s with other characterwere found to have very little impact, thus these results are not shown here. Time [a.u.] −50−2502550 C u rr e n d e n s i y [ − a . u . ] index of lo se 12345 678910 0 2 4 6 8 10 Index of lo se −5 −4 R M S E [ a . u . ] Figure 7.
Convergence of the current density in silicon with increasing number of lo’s(left). The following lo settings have been considered (number of lo’s corresponding toangular momenta in parentheses) 1 (2 s, 1 p); 2 (2 s, 2 p) 3 (2 s, 2 p, 1 d); 4 (2 s, 2 p,2 d); 5 (3 s, 2 p, 2 d); 6 (3 s, 3 p, 2 d); 7 (3 s, 3 p, 3 d); 8 (3 s, 4 p, 3 d); 9 (4 s, 4 p, 3d); 10 (4 s, 4 p, 4 d). Right: corresponding RMSE.
3. Benchmark results
In this section, we present a benchmark of our implementation, contrasting theimaginary part of the dielectric function obtained with Eq. (10) with that of theLR-TDDFT, employing the adiabatic local-density approximation (ALDA) as alreadyimplemented in exciting [25]. We also compare the current density obtained with ourimplementation with results from Octopus [26].
As prototypical materials for our initial benchmark, we choose silicon and 2-dimensionalMoS . The imaginary part of their dielectric functions are given in Figs. 8 and 9,respectively. In the RT-TDDFT calculations, we considered an impulsive electric fieldwith an amplitude small enough to not induce deviations from the linear regime. Apartfrom the already commented “fake” plasmon in RT-TDDFT for smaller frequencies, weobserve overall a remarkable agreement between both results. A similar comparison fordiamond is provided in Fig. C7. Octopus has been one of the first codes to evaluate the propagation of KS wavefunctionswithin the framework of RT-TDDFT [3, 27, 32, 26]. Hence, it is a most suitable package ll-electron full-potential implementation of real-time TDDFT in exciting Figure 8.
Imaginary part of the dielectric function of bulk silicon: Comparisonbetween RT- and LR-TDDFT.
Figure 9.
Out-of-plane ( zz ) (left) and in-plane ( xx ) (right) tensor-components ofthe imaginary part of the dielectric function of 2-dimensional MoS obtained by RT-TDDFT in comparison with LR-TDDFT results .to benchmark our results, even though it employs a different scheme to solve the KSequations, namely pseudo-potentials combined with a real-space mesh [26]. On the leftside of Fig. 10, we depict the current density in cubic BN as response to an impulsiveelectric field. The agreement between the results of exciting and Octopus is impressive.On the right side, we provide the imaginary part of the dielectric function. This servesas well as a measure of how similar the Fourier-transforms of both curves are. We showalso the result from our LR-TDDFT calculation, depicted as gray-shaded area. Onceagain, the agreement is excellent.As a second benchmark, we depict in Fig. 11 the current density in Si exposedto an electric field, whose expression follows Eq. (17) and is shown in the inset, andcompare the result to that obtained with Octopus. An interesting aspect here is thatthis field is strong enough to induce a nonlinear response, as it can be seen from theresidual current density (after t = 700 a.u., when the external field turns to zero). Also ll-electron full-potential implementation of real-time TDDFT in exciting Figure 10.
Comparison between the current density obtained with exciting andOctopus for an impulsive electric field applied to cubic BN. The right panel comparescorresponding results for the imaginary part of dielectric function to that obtainedwith the LR-TDDFT implementation of exciting (gray shaded area). in this nonlinear regime, the agreement between exciting and Octopus is very good.Similar agreement is found for SiC, see Fig. C8.
Time [a.u.] −0.00010.00000.00010.00020.0003 C u rr e n t d e n s i t y [ a . u . ] exciting Octopus
Time [a.u.] −2.50.02.5 𝘌 ( 𝘵 ) Figure 11.
Comparison between the current density in Si, obtained with exciting and Octopus. The inset depicts the applied electric field (in 10 a.u.).
4. Implemented features
The RT-TDDFT implementation naturally provides the evolution of the KS system, i.e. , KS energies and wavefunctions, charge density, and total energy as functions oftime. In this section, we choose three features of our implementation which highlight itas an interesting tool to aid the interpretation of experiments. ll-electron full-potential implementation of real-time TDDFT in exciting We start with the dynamics of an excitation in two-dimensional MoS caused by a laserpulse with the electric field given by Eq. (17) and plotted in the inset on the rightpanel of Fig. 12. The pulse duration is set to T pulse = 400 a.u., the frequency to ω = 0 .
15 a.u. (corresponding to a photon energy of 4 .
08 eV), and the peak intensityto E m = 0 . x direction), the other one perpendicular to it ( z direction).Figure 12 shows the current density in these two cases on the left, as well as the numberof excited electrons on the right. Time [a.u.] −600−400−2000200400600 C u rr e n t d e n s i t y [ − a . u . ] x z 0 100 200 300 400 Time [a.u.] 𝘕 𝘦 𝘹 𝘤 ( 𝘵 ) 𝘵 [a.u.] −0.010.000.01 𝘌 ( 𝘵 ) [ a . u . ] Figure 12.
Left: Current density in MoS under the action of electric fields paralleland perpendicular to the monolayer plane ( x and z directions, respectively). Right:number of electrons per unit cell excited to the conduction band. The inset shows thetime dependence of the applied field. The current density is considerably higher in the case of in-plane polarization witha peak height being about 3 times larger. Some nonlinear effects are already observable.When we compare the number of excited electrons, we observe that, in the end, after thepulse is removed, 1.26 electrons per unit cell remain excited in the case of polarizationalong the x direction, but two orders of magnitude less, i.e. , 0.012, for the z direction.We can understand this difference by the 2D nature of the material and can trace itback to the dielectric function (Fig. 9). The out-of-plane component ε zz at 4 .
08 eV ismuch higher than the in-plane component which means that, at this frequency, MoS can absorb electromagnetic waves with the electric field parallel to the monolayer planemuch better than perpendicular to it.It is also interesting to observe from and to which bands the electrons are excited.In Fig. 13 we show for three different times, i.e. , t = 100, 200 and 400 a.u. how theexcitations are distributed over the k -space. In the top panels, we provide the resultsfor the x -polarization and in the bottom panels for the z -polarization. In the case ofin-plane polarization, some excitations are present at t =100 a.u.; many more electronsbecome excited at t =200 a.u., followed by a decrease thereafter. Interestingly, around ll-electron full-potential implementation of real-time TDDFT in exciting M K−505 E n e r g y [ e V ] 𝘵 = 100 a.u. M K𝘵 = 200 a.u. M K𝘵 = 400 a.u.M K−505 E n e r g y [ e V ] M K M K
Figure 13.
Band structure of MoS along the ΓM and ΓK directions as the responseto a laser-pulse with the electric field along (top panels) and perpendicular to (bottompanels) the MoS -plane. The circles in red (green) indicate the degree of population(depopulation) of the conduction (valence) bands at the specified times. the Γ and K points, the holes tend to be formed not on the valence-band top, but indeeper-lying bands, whereas at the M point, the holes are predominantly at the topone. In the case of the perpendicular polarization, there is almost no difference betweenthe excitations at times t =100 a.u. and t =200 a.u. In the end, only a few excitationsremain, and they are not concentrated at the band edges, but rather in deeper- andhigher-lying bands, respectively. We now analyze the response of silicon exposed to an electric field along the [001]direction whose expression follows Eq. (17), with the parameters ω = 0 . T pulse = 744 a.u. = 18.0 fs, and t = 0,according to Ref. [9]. The amplitude E m is varied so that we can observe a progressionfrom the linear to the nonlinear regime.In Fig. 14, we depict the Fourier-transform of the polarization field P =( D − E ) / (4 π ) normalized by the amplitude E m of the applied electric field. Theintensity of the electromagnetic wave I (in W/cm ) is obtained from the amplitudeas I = 3 . × E m . It is important to recall that E m stands for the amplitudeof the electric field inside the bulk material. Due to reflection at the interface andboundary conditions, the electric field generated by the exciting laser may be differentfrom E m , sometimes even two orders of magnitude higher [9]. We can observe that,for intensities of 9 . × W/cm and higher, the third and even the fifth harmoniccomponents are excited, and these components are obviously stronger the more intensethe field is.We also evaluate the number of excited electrons per unit cell, n ex , at a sufficient ll-electron full-potential implementation of real-time TDDFT in exciting / −5 −3 −1 𝘗 ( ) / 𝘌 𝘮 Intensity [10 W/cm ]0.061 624 96217 386603 Figure 14.
Fourier transform of the polarization field in silicon under the actionof a laser pulse with fundamental frequency ω = 0 . Intensity [W/cm²] −6 −4 −2 N u m b e r o f e x c i t e d e l e c t r o n s Figure 15.
Number of excited electrons per unit cell in silicon as a function of theintensity of an external electric pulse of frequency ω = 0 . large time after the electric field has been switched off, as shown in Fig. 15. We recognizethat n ex is connected to the intensity I of the electromagnetic wave by a power law, i.e. , n ex = CI n . (18)By means of a least square fit, we find n = 1 .
94, which agrees with Ref. [42].
We now simulate a pump-probe experiment, taking diamond as test material. Theelectric field of the pump pulse, given by Eq. (17), has a gaussian-like envelope withwidth T pulse = 644 a.u., amplitude E m = 2 . ω =0 . t = 700 a.u., a weak impulsive electric field E ( t ) = 0 . δ ( t − J pump − probe ) and the pump pulse J pump , as obtained from two separate calculations. ll-electron full-potential implementation of real-time TDDFT in exciting Energy [eV] I m ( ) WithpumpWithoutpump
Figure 16.
Imaginary part of the dielectric function of diamond probed after anelectric field acting had been applied as pump (blue curve). For comparison, the curveexpected without a pumping field is given in green.
5. Conclusions
In this paper, we have presented the implementation of RT-TDDFT in the full-potentialLAPW+lo package exciting , providing the underlying theory and details on theconvergence behavior as well as parallelization performance. As benchmarks, we havecompared our results with those obtained with the Octopus code as well as LR-TDDFTresults from exciting , finding excellent agreement in all cases. We have shown threeexamples of applications how our implementation could be used for the interpretationof experiments. These are the excitation dynamics of a material upon radiation witha laser pulse, the non-linear response of a material to laser pulses, and the dielectricfunction after a pump pulse. The implementation is included in the latest release, exciting oxygen, and the code can be downloaded for free from the exciting webpage[24]. All data presented here are available in the NOMAD Repository [44, 45] (DOI:10.17172/NOMAD/2021.01.20-1).
Acknowledgments
This work was supported by the Deutsche Forschungsgemeinschaft (DFG)- Projektnum-mer 182087777 - SFB 951. We thank Keith Gilmore, Santiago Rigamonti, Sven Lubeckand, Felix Henneke for the critical review of this manuscript. Alexander Buccheri andSebastian Tillack are acknowledged for reviewing our code. ll-electron full-potential implementation of real-time TDDFT in exciting References [1] Runge E and Gross E K U 1984
Phys. Rev. Lett. (12) 997–1000[2] Botti S, Schindlmayr A, Sole R D and Reining L 2007 Reports on Progress in Physics
357 ISSN0034-4885[3] Andrade X, Alberdi-Rodriguez J, Strubbe D A, Oliveira M J T, Nogueira F, Castro A, MuguerzaJ, Arruabarrena A, Louie S G, Aspuru-Guzik A, Rubio A and Marques M A L 2012
Journal ofPhysics: Condensed Matter Phys. Chem. Chem. Phys. (40) 26599–26606[5] Maitra N T 2016 The Journal of Chemical Physics
Journal of Chemical Theory andComputation The Journal of Chemical Physics
Annual Reports in Computational Chemistry ( AnnualReports in Computational Chemistry vol 11) ed Dixon D A (Elsevier) pp 103 – 146[9] Yabana K, Sugiyama T, Shinohara Y, Otobe T and Bertsch G F 2012
Phys. Rev. B (4) 045134[10] Meng S and Kaxiras E 2008 The Journal of Chemical Physics
The Journal of Chemical Physics
Advanced Theory and Simulations The Journal of Chemical Physics
Physical Review B Applied Sciences Preprint )[16] Miyamoto Y, Zhang H, Miyazaki T and Rubio A 2015
Physical Review Letters
Journal of Physics: Condensed Matter Phys. Rev. B (8)081204[19] Pemmaraju C D 2019 Computational Condensed Matter e00348 ISSN 2352-2143[20] Imamura Y, Suzuki K, Iizuka T and Nakai H 2015 Chemical Physics Letters
30 – 36 ISSN0009-2614[21] Rigamonti S, Botti S, Veniard V, Draxl C, Reining L and Sottile F 2015
Phys. Rev. Lett. (14)146402[22] Lejaeghere K, Bihlmayer G, Bj¨orkman T, Blaha P, Bl¨ugel S, Blum V, Caliste D, Castelli I E,Clark S J, Dal Corso A, de Gironcoli S, Deutsch T, Dewhurst J K, Di Marco I, Draxl C,Du lak M, Eriksson O, Flores-Livas J A, Garrity K F, Genovese L, Giannozzi P, GiantomassiM, Goedecker S, Gonze X, Gr˚an¨as O, Gross E K U, Gulans A, Gygi F, Hamann D R, HasnipP J, Holzwarth N A W, Iu¸san D, Jochym D B, Jollet F, Jones D, Kresse G, Koepernik K,K¨u¸c¨ukbenli E, Kvashnin Y O, Locht I L M, Lubeck S, Marsman M, Marzari N, Nitzsche U,Nordstr¨om L, Ozaki T, Paulatto L, Pickard C J, Poelmans W, Probert M I J, Refson K, RichterM, Rignanese G M, Saha S, Scheffler M, Schlipf M, Schwarz K, Sharma S, Tavazza F, Thunstr¨omP, Tkatchenko A, Torrent M, Vanderbilt D, van Setten M J, Van Speybroeck V, Wills J M, YatesJ R, Zhang G X and Cottenier S 2016
Science
ISSN 0036-8075[23] Gulans A, Kozhevnikov A and Draxl C 2018
Phys. Rev. B (16) 161105[24] The exciting code http://exciting-code.org/ ll-electron full-potential implementation of real-time TDDFT in exciting [25] Gulans A, Kontur S, Meisenbichler C, Nabok D, Pavone P, Rigamonti S, Sagmeister S, Werner Uand Draxl C 2014 Journal of Physics: Condensed Matter The Journal of Chemical Physics
Phys. Rev. B (12) 7998–8002[28] Yabana K, Nakatsukasa T, Iwata J and Bertsch G F 2006 Phys. Stat. Sol. (b)
Computer PhysicsCommunications
30 – 38 ISSN 0010-4655[30] Marques M A L, Ullrich C A, Nogueira F, Rubio A, Burke K and Gross E K U (eds) 2006
Time-Dependent Density Functional Theory
Lecture Notes in Physics, 706 (Berlin, Heidelberg:Springer) ISBN 9783540354260[31] Marques M A L, Maitra N T, Nogueira F M S, Gross E K U and Rubio A (eds) 2012
Fundamentals of Time-Dependent Density Functional Theory
Lecture Notes in Physics, 837(Berlin, Heidelberg: Springer) ISBN 9783642235184[32] Castro A, Marques M A L and Rubio A 2004
The Journal of Chemical Physics
Journal of Chemical Theory andComputation Computer Physics Communications
92– 95 ISSN 0010-4655[35] Vorwerk C, Aurich B, Cocchi C and Draxl C 2019
Electronic Structure Introduction to Electrodynamics (Harlow: Pearson Education UK) ISBN129202142X[37] Temnov V V, Sokolowski-Tinten K, Zhou P, El-Khamhawy A and von der Linde D 2006
Phys.Rev. Lett. (23) 237403[38] Sato S A, Yabana K, Shinohara Y, Otobe T and Bertsch G F 2014 Physical Review B Phys. Rev. B (4) 2643–2650[40] Zheng Q, Chu W, Zhao C, Zhang L, Guo H, Wang Y, Jiang X and Zhao J 2019 WileyInterdisciplinary Reviews: Computational Molecular Science e1411 ISSN 1759-0876[41] Li Y, He S, Russakoff A and Varga K 2016
Physical Review E Phys. Rev. B (16) 165104[43] Otobe T, Yabana K and Iwata J I 2009 Journal of Physics: Condensed Matter Journal of Physics: Materials https://doi.org/10.1088/2515-7639/ab13bb [45] Nomad repository, dataset: rttddft-exciting. https://dx.doi.org/10.17172/NOMAD/2021.01.20-1/ Appendix A. Input
In Fig. A1, we display as an example, the input file of MoS .For the RT-TDDFTcalculations, the most important elements are captured by the element rt tddft within the excited-state module xs . Here, the input file defines as propagator CFM4(commutator free Magnus of 4th order), with an evolution time starting at 0 (default)up to 400 a.u. ( endtime ), with steps of 0 .
05 a.u. ( timestep ). The vector potential ll-electron full-potential implementation of real-time TDDFT in exciting Figure A1.
Example of input file ( input.xml ). A ( t ), assuming the transverse geometry ( afield="total" ), is described by the element laser . In this case, we have a field applied along the z axis, with gaussian-like envelopeas in Eq. (17), where A m =100 a.u., ω = 0 .
15 a.u., T pulse = 400 a.u., t = 0 a.u., anda null extra phase for the term cos( ω t ). The element symmetry break defines an axis,given in Cartesian coordinates ( cartesian="true" ), to break the crystal symmetry.The full reference of input variables and their meaning is provided at the exciting webpage [24]. ll-electron full-potential implementation of real-time TDDFT in exciting Appendix B. Propagators
In this section, expressions for the propagator ˆ U ( t + ∆ t, t ) are provided. The derivationsand assessments of their advantages or disadvantages can be found in Refs. [32], [33],and [34]. We start with the most basic extension to the simple exponential propagator,namely the exponential at the midpoint:ˆ U ( t + ∆ t, t ) = exp (cid:20) − i∆ t ˆ H (cid:18) t + ∆ t (cid:19)(cid:21) , (B.1)where, the extrapolation for obtaining the Hamiltonian ˆ H ( t + f ∆ t ) at future times is:ˆ H ( t + f ∆ t ) = (1 + f ) ˆ H ( t ) − f ˆ H ( t − ∆ t ) . (B.2)Another extension to the simple exponential propagator which keeps time-reversalsymmetry to be approximately fulfilled isˆ U ( t + ∆ t, t ) = exp (cid:20) − i ∆ t H ( t + ∆ t ) (cid:21) exp (cid:20) − i ∆ t H ( t ) (cid:21) . (B.3)Further improvement is provided by the so called commutator-Free Magnus expansionof 4th order: ˆ U ( t + ∆ t, t ) = exp h − i∆ t ( α ˆ H ( t ) + α ˆ H ( t )) i ×× exp h − i∆ t ( α ˆ H ( t ) + α ˆ H ( t )) i , (B.4)where t , = t + (cid:16) ∓ √ (cid:17) ∆ t , and α , = ∓ √ .A different approach is evaluating the exponential operator exactly rather thanTaylor-expanding it. This can be done by taking into account an adiabatic basis formedby the eigenvectors of ˆ H ( t ). It means that for each t , we solveˆ H ( t ) | φ j k ( t ) i = ε j k ( t ) | φ j k ( t ) i (B.5)and then expand | ψ j k ( t ) i = X i α ij k ( t ) | φ i k ( t ) i , (B.6)where α ij k ( t ) = h φ k i ( t ) | ψ k j ( t ) i . Sinceˆ U ( t + ∆ t, t ) | φ m k ( t ) i = e − i ε m k ( t )∆ t | φ m k ( t ) i , (B.7)when considering ˆ U ( t + ∆ t, t ) in the form of the simple exponential propagator, theaction of the propagator, using Eq. (B.6), isˆ U ( t + ∆ t, t ) | ψ j k ( t ) i = X i α ij k ( t )e − i ε m k ( t )∆ t | φ m k ( t ) i . (B.8)We further utilize then the expansion of | φ m k ( t ) i in terms of our LAPW+lo basis.Although the exponential operator is exactly obtained, i.e. , without the need of aTaylor expansion, this approach now relies on an expansion in terms of the adiabaticbasis and on the assumption of a simple exponential for the propagator. A first ll-electron full-potential implementation of real-time TDDFT in exciting H ( t + ∆ t/ i.e. , | ψ j k ( t + ∆ t ) i = | ψ j k ( t ) i − i∆ t k + 2 k + 2 k + k ) (B.9)where k = S − k H k ( t ) | ψ j k ( t ) i , (B.10) k = S − k H k (cid:18) t + ∆ t (cid:19) (cid:20) | ψ j k ( t ) i + k ∆ t (cid:21) , (B.11) k = S − k H k (cid:18) t + ∆ t (cid:19) (cid:20) | ψ j k ( t ) i + k ∆ t (cid:21) , (B.12) k = S − k H k ( t + ∆ t ) [ | ψ j k ( t ) i + k ∆ t ] . (B.13) Appendix C. Additional results
Appendix C.1. Convergence behaviorAppendix C.1.1. Time step
In Figs. C1 and C2, we show the convergence behaviorof the current density in silicon, probing the size of the time step after an external
Time [a.u.] −0.0010.0000.0010.0020.003 C u rr e n t d e n s i t y [ a . u . ] Time step [a.u.]0.0010.0020.0050.010.020.050.10 10 −3 −2 −1 Time step [a.u.] −4 −3 R M S E [ a . u . ] 𝘵 𝘮 − 𝘵 𝘵 𝘮 𝘵 𝘮 + 𝘵𝘑 𝘮 − 𝘑𝘑 𝘮 𝘑 𝘮 + 𝘑 Figure C1.
Convergence behavior of the current density in silicon with respectto time-step when an impulsive external field is applied (left). The inset amplifiesthe region around the first maximum t m =41.5 a.u., with J m = 9 . × − a.u., δJ = 1 × − a.u., and δt = 0 . electric field is applied along the [001] direction. In the first case, it is a delta function, E ( t ) = 0 . δ ( t − .
5) a.u., while in the second case, it has a gaussian-like envelop, Eq.(17), with the parameters E m = 4 . × − , ω = 0 . t = 2 .
0, and T pulse = 452 (inatomic units). ll-electron full-potential implementation of real-time TDDFT in exciting Time [a.u.] −50−2502550 C u rr e n t d e n s i t y [ − a . u . ] Time Step [a.u.]0.010.020.050.100.20 10 −2 −1 Time step [a.u.] −7 −6 R M S E [ a . u . ] Figure C2.
Convergence behavior of the current density in silicon with respect totime-step when an external field with gaussian-like envelope is applied. The RMSEis displayed on the right, taking the results obtained with a time-step 0 .
01 a.u.asreference.
Appendix C.1.2. k -points Figures C3 and C4 display the convergence behavior of thecurrent density in silicon with respect to the k -grid. The external electric fields are thesame as in Appendix C.1.1. Time [a.u.] −0.0010.0000.0010.0020.003 C u rr e n d e n s i y [ a . u . ] Number of 𝗸-points24 Number of 𝗸-points −5 −4 −3 R M S E [ a . u . ] Figure C3.
Impact of the number of k -points on the current density in silicon, whenan impulsive electric field is applied. The right panel shows the RMSE, taking thecalculation with 24 × ×
24 as reference.
Figure C5 shows the influence of the number of k -points on the dielectric functionof carbon. These results have been obtained from a Fourier transform of those given inFig. 3 (calculations up to a time of 5000 a.u.). ll-electron full-potential implementation of real-time TDDFT in exciting Time [a.u.] −50050 C u rr e n d e n s i y [ − a . u . ] Number of 𝗸-points20 Number of 𝗸-points −6 −5 R M S E [ a . u . ] Figure C4.
Convergence behavior of the current density in silicon with respect to the k -points, where the external field has a gaussian-like envelope. The right panel showsthe RMSE taking the calculation with highest a k -grid of 20 × ×
20 as reference.
Energy [eV] I m ( ) Number of 𝗸-points14 Energy [eV] −10−5051015 R e ( ) Figure C5.
Impact of the number of k -points on the convergence of the imaginary(left) and real (right) parts of the dielectric function of diamond. Appendix C.1.3. Basis-set size
Figure C6 displays the impact of the choice of theparameter rgkmax on the convergence behavior of the current density in silicon. Theexternal electric field has a gaussian-like envelope, as given by Eq. (17), with sameparameters as in Appendix C.1.1.
Appendix C.2. Benchmarks complementing Section 3
Figure C7 displays how the imaginary part of the dielectric function of diamond obtainedwith RT-TDDFT compares with that from LR-TDDFT. In Fig. C8, we compare thecurrent density in SiC obtained with exciting to the result by Octopus. The externalfield is a delta function applied along the [001] direction, i.e., D ( t ) = 0 . δ ( t −
2) [a.u.].The right side of the figure shows how the imaginary part of the dielectric function ll-electron full-potential implementation of real-time TDDFT in exciting Time [a.u.]
60 40 200204060 C u rr e n t d e n s i t y [ − a . u . ] rgkmax7.06.56.05.55.04.54.0 4 5 6 7 rgkmax [a.u.] R M S E [ a . u . ] Figure C6.
Convergence behavior of the current density in silicon with respect to thebasis-set size, determined by the parameter rgkmax . The external electric field has agaussian-like envelop. The right panel depicts the RMSE, taking the calculation withhighest rgkmax as reference.
Figure C7.
Imaginary part of the dielectric function of diamond: Comparisonbetween RT- and LR-TDDFT. obtained from the RT-TDDFT calculations of the two codes compares with LR-TDDFTobtained with exciting . ll-electron full-potential implementation of real-time TDDFT in exciting Figure C8.
Current density in SiC, exposed to an impulsive field, as obtained with exciting compared to the results of Octopus. The right panel compares the imaginarypart of the dielectric function obtained from the RT-TDDFT implementations of bothcodes to the LR-TDDFT result obtained with excitingexciting