Accurate time propagation method for the coupled Maxwell and Kohn-Sham equations

An accurate method for time propagation of the coupled Maxwell and time-dependent Kohn-Sham (TDKS) equation is presented. The new approach uses a simultaneous fourth-order Runge-Kutta–based propagation of the vector potential and the Kohn-Sham orbitals. The approach is compared to the conventional fourth-order Taylor propagation and predictor-corrector methods. The calculations show several computational and numerical advantages, including higher computational performance, greater stability, better accuracy, and faster convergence.

Figure 1

(a) Current and (b) vector potential in the $z$ direction in a diamond crystal induced by a $\delta$ kick applied in the $z$ direction. A $2\times 2\times 2 k$-point mesh is used. The Kohn-Sham orbitals are propagated with the Taylor propagation scheme (Algorithm 1) for times up to 400 a.u. Current plotted only up to 200 a.u. Well converged results are obtained for time steps $\mathrm{\Delta }t\le 0.002$ a.u.

Figure 2

(a) Current and (b) vector potential in the $z$ direction in a diamond crystal induced by a $\delta$ kick applied in the $z$ direction. A $2\times 2\times 2 k$-point mesh is used. The Kohn-Sham orbitals are propagated with the predictor-corrector propagation scheme (Algorithm 2) for times up to 400 a.u. Current plotted only up to 200 a.u. The propagation only remains stable for very small time steps, $\mathrm{\Delta }t\le 0.002$ a.u., with the Taylor benchmark and PC method exactly overlapping.

Figure 3

(a) Ground state density isosurface plot and (b) number of electrons in the ground state orbitals as a function of time.

Figure 4

(a) Current, (b) vector potential, and (c) total energy change of a diamond crystal after an applied delta kick on a $2\times 2\times 2 k$-point mesh. The system is propagated up to a time of 400 a.u. The current is only plotted up to 200 a.u. The PC and SRK4 propagation schemes are compared. $\mathrm{\Delta }{t}_{PC}$ = 0.02 a.u. is used for the PC method, in agreement with the previous sections. The results of the benchmark Taylor and PC methods, also shown in Fig. 2, are shown here for comparison. $\mathrm{\Delta }{t}_{\text{SRK}4}$ = 0.05 a.u. is used since it is the maximum allowed value for this scheme. The SRK4 method finds excellent agreement with the benchmark Taylor propagation with small time step, $\mathrm{\Delta }{t}_T$ = 0.002 a.u.

Figure 5

(a) Current and (b) vector potential of a diamond crystal after an applied $\delta$ kick on a $5\times 5\times 5 k$-point mesh. The system is propagated up to a time of 2000 a.u. The current is plotted up from 1700 to 2000 a.u. Inset (a) shows the current in units of ${10}^{-3}$ a.u. from time 1050 to 1150 a.u. Inset (b) shows the $z$ vector potential in a.u. from time 100 to 250 a.u. The PC and SRK4 propagation schemes are compared. Time steps of 0.02 and 0.05 a.u. are shown for the PC method. A time step of 0.05 a.u. is used for the SRK4 method. The SRK4 method finds excellent agreement with the benchmark Taylor propagation with small time step, $\mathrm{\Delta }{t}_T$ = 0.002 a.u. for the duration of the propagation.

Figure 6

Inverse of the dielectric constant of a diamond crystal obtained through TDDFT simulations of the electron dynamics after a delta kick on a $5\times 5\times 5 k$-point mesh. The Kohn-Sham orbitals are propagated with the Taylor, PC, and SRK4 propagation schemes. Plots (a) and (b) show the (a) real and (b) imaginary parts of the inverse of the dielectric constant obtained by Fourier transforming the induced vector potential with a small broadening constant, $\eta =0.005$ a.u. Plots (c) and (d) show the (c) real and (d) imaginary parts of the inverse of the dielectric constant obtained without a broadening parameter.

Figure 7

Real and imaginary components of the dielectric function of graphene obtained with TDDFT simulations of the electron dynamics after a $\delta$ kick on a $11\times 19\times 1 k$-point mesh. The Kohn-Sham orbitals are propagated with the SRK4 scheme with time step, 0.05 a.u.

Figure 8

(a) Electric field and current in diamond subject to a short laser pulse, calculated by the SRK4 method on a $12\times 12\times 12 k$-point mesh. The pulse is 1240 a.u. (30 fs) wide with 0.057 a.u. (1.55 eV) frequency and 0.0154 a.u. amplitude. (b) The first three harmonic generators located at 0.057, 0.171, and 0.285 a.u., as shown in the logarithmic scaled current transformed in energy space.