^{1}, Erik W. Draeger

^{2}, Yosuke Kanai

^{1,3,a)}and Alfredo A. Correa

^{1,a)}

### Abstract

Explicit integrators for real-time propagation of time-dependent Kohn-Sham equations are compared regarding their suitability for performing large-scale simulations. Four algorithms are implemented and assessed for both stability and accuracy within a plane-wave pseudopotential framework, employing the adiabatic approximation to the exchange-correlation functional. Simulation results for a single sodium atom and a sodium atom embedded in bulk magnesium oxide are discussed. While the first-order Euler scheme and the second-order finite-difference scheme are unstable, the fourth-order Runge-Kutta scheme is found to be conditionally stable and accurate within this framework. Excellent parallel scalability of the algorithm up to more than a thousand processors is demonstrated for a system containing hundreds of electrons, evidencing the suitability for large-scale simulations based on real-time propagation of time-dependent Kohn-Sham equations.

We gratefully acknowledge fruitful discussions with E. K. U. Gross, A. Rubio, K. Burke, N. Maitra, and J. Kohanoff. A.S. is supported through the Physical and Life Sciences Directorate at Lawrence Livermore National Laboratory. Computing support for this work came from the Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge program. Part of this work was performed under the auspices of the U.S. Department of Energy at Lawrence Livermore National Laboratory under Contract No. DE-AC52-07A27344.

I. INTRODUCTION

II. THEORETICAL FRAMEWORK AND NUMERICAL TREATMENT

A. Time-dependent Kohn-Sham equations

B. Energy conservation during the propagation

C. Plane-wave expansion of the wave functions

D. Computational details

III. PROPAGATION SCHEMES

A. General considerations

B. First-order Euler scheme

C. Second-order finite differences

D. Second-order Runge-Kutta scheme

E. Fourth-order Runge-Kutta scheme

IV. TEST CASE I: THE ISOLATED NA ATOM

A. Numerical stability and conservation of energy

B. Evolution of the system

V. TEST CASE II: NA ATOM EMBEDDED IN MGO

A. Parallel scaling

VI. CONCLUSIONS

### Key Topics

- Wave functions
- 31.0
- Sodium
- 12.0
- Ground states
- 9.0
- Orthogonalization
- 6.0
- Hohenberg Kohn theorem
- 5.0

## Figures

Cubic unit cell containing a single Na atom (small circle). To prepare a non-equilibrium initial condition, the Na 3*s* wave function (represented by the yellow isosurface) has been shifted by (0.32, 0.32, 0.32) Å, see (a). In (b) and (c), two snapshots of the RK4 simulation are shown.

Cubic unit cell containing a single Na atom (small circle). To prepare a non-equilibrium initial condition, the Na 3*s* wave function (represented by the yellow isosurface) has been shifted by (0.32, 0.32, 0.32) Å, see (a). In (b) and (c), two snapshots of the RK4 simulation are shown.

Total energy *E* _{tot} (in eV) of the Na atom (at *t* = 0 fs, the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). In (a), the Euler scheme (black solid line, Δ*t* = 0.069 as) and the sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) are compared to the Runge-Kutta propagators (blue solid line). The second-order (Δ*t* = 0.069 as) and the fourth-order (Δ*t* = 0.691 as) Runge-Kutta scheme yield the same trajectory for the times shown in (a). In (b), the fully self-consistent second-order finite-difference method (green solid line) is compared to the sc100-SOD (red solid line) and the non-self-consistent (red dotted line) one for Δ*t* = 0.069 as.

Total energy *E* _{tot} (in eV) of the Na atom (at *t* = 0 fs, the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). In (a), the Euler scheme (black solid line, Δ*t* = 0.069 as) and the sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) are compared to the Runge-Kutta propagators (blue solid line). The second-order (Δ*t* = 0.069 as) and the fourth-order (Δ*t* = 0.691 as) Runge-Kutta scheme yield the same trajectory for the times shown in (a). In (b), the fully self-consistent second-order finite-difference method (green solid line) is compared to the sc100-SOD (red solid line) and the non-self-consistent (red dotted line) one for Δ*t* = 0.069 as.

Total energy *E* _{tot} (in eV) of the Na atom as a function of time *t* (in fs). The curves in (a) result from the sc100-SOD second-order finite-difference scheme and the ones in (b) from the second-order Runge-Kutta scheme. Time steps of Δ*t* = 0.069 as (black solid lines), Δ*t* = 0.104 as (red solid lines), and Δ*t* = 0.138 as (green solid lines) were used.

Total energy *E* _{tot} (in eV) of the Na atom as a function of time *t* (in fs). The curves in (a) result from the sc100-SOD second-order finite-difference scheme and the ones in (b) from the second-order Runge-Kutta scheme. Time steps of Δ*t* = 0.069 as (black solid lines), Δ*t* = 0.104 as (red solid lines), and Δ*t* = 0.138 as (green solid lines) were used.

Time steps Δ*t* (in as) for which the non-self-consistent second-order finite-difference scheme (a) or the fourth-order Runge-Kutta scheme (b) are stable (green triangles pointing up) or unstable (red triangles pointing down), depending on the respective plane-wave cutoff energy *E* _{cut} (in Ry) used.

Time steps Δ*t* (in as) for which the non-self-consistent second-order finite-difference scheme (a) or the fourth-order Runge-Kutta scheme (b) are stable (green triangles pointing up) or unstable (red triangles pointing down), depending on the respective plane-wave cutoff energy *E* _{cut} (in Ry) used.

Sum of the expectation values (in eV) of all valence wave functions *i* of the Na atom (at *t* = 0 fs the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). In (a) the Euler scheme (black solid line, Δ*t* = 0.069 as) and the sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) are compared to the Runge-Kutta propagators (blue solid line). The second-order (Δ*t* = 0.069 as) and the fourth-order (Δ*t* = 0.691 as) Runge-Kutta scheme yield the same trajectory for the times shown in (a). In (b) the fully self-consistent second-order finite-difference method (green solid line) is compared to the sc100-SOD (red solid line) and the non-self-consistent (red dotted line) one for Δ*t* = 0.069 as.

Sum of the expectation values (in eV) of all valence wave functions *i* of the Na atom (at *t* = 0 fs the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). In (a) the Euler scheme (black solid line, Δ*t* = 0.069 as) and the sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) are compared to the Runge-Kutta propagators (blue solid line). The second-order (Δ*t* = 0.069 as) and the fourth-order (Δ*t* = 0.691 as) Runge-Kutta scheme yield the same trajectory for the times shown in (a). In (b) the fully self-consistent second-order finite-difference method (green solid line) is compared to the sc100-SOD (red solid line) and the non-self-consistent (red dotted line) one for Δ*t* = 0.069 as.

Total energy *E* _{tot} (in eV) of the Na atom (at *t* = 0 fs the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). The results have been obtained using the RK4 propagator and Δ*t* = 0.691 as.

Total energy *E* _{tot} (in eV) of the Na atom (at *t* = 0 fs the 3*s* wave function was shifted by (0.32, 0.32, 0.32) Å from its equilibrium position in real space) as a function of time *t* (in fs). The results have been obtained using the RK4 propagator and Δ*t* = 0.691 as.

The 64-atom unit cell of MgO (Mg atoms red circles, O atoms blue circles) containing a single Na atom (black circle) on an oxygen lattice position.

The 64-atom unit cell of MgO (Mg atoms red circles, O atoms blue circles) containing a single Na atom (black circle) on an oxygen lattice position.

Sum of the expectation values (in eV) of all valence wave functions in the Na:MgO 64-atom supercell (at *t* = 0 fs the Na-induced level within the MgO gap was shifted by (0.032, 0.032, 0.032) Å from its equilibrium position in real space) as a function of time *t* (in fs). The sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) is compared to the second-order (green solid line, Δ*t* = 0.069 as) and the fourth-order (blue solid line, Δ*t* = 0.691 as) Runge-Kutta scheme.

Sum of the expectation values (in eV) of all valence wave functions in the Na:MgO 64-atom supercell (at *t* = 0 fs the Na-induced level within the MgO gap was shifted by (0.032, 0.032, 0.032) Å from its equilibrium position in real space) as a function of time *t* (in fs). The sc100-SOD second-order finite-difference scheme (red solid line, Δ*t* = 0.069 as) is compared to the second-order (green solid line, Δ*t* = 0.069 as) and the fourth-order (blue solid line, Δ*t* = 0.691 as) Runge-Kutta scheme.

Total energy *E* _{tot} (in eV) of the Na:MgO 64-atom supercell (at *t* = 0 fs the Na-induced level within the MgO gap was shifted by (0.032, 0.032, 0.032) Å from its equilibrium position in real space) as a function of time *t* (in fs). The results have been obtained using the fourth-order Runge-Kutta propagator and Δ*t* = 0.691 as.

Total energy *E* _{tot} (in eV) of the Na:MgO 64-atom supercell (at *t* = 0 fs the Na-induced level within the MgO gap was shifted by (0.032, 0.032, 0.032) Å from its equilibrium position in real space) as a function of time *t* (in fs). The results have been obtained using the fourth-order Runge-Kutta propagator and Δ*t* = 0.691 as.

The number of steps that can be performed within 1 s wall time is plotted versus the number of processing cores used in the calculation for the Na:MgO 64-atom supercell. The steepest descent algorithm (black curve) is compared to the RK4 propagation (red curve).

The number of steps that can be performed within 1 s wall time is plotted versus the number of processing cores used in the calculation for the Na:MgO 64-atom supercell. The steepest descent algorithm (black curve) is compared to the RK4 propagation (red curve).

Article metrics loading...

Full text loading...

Commenting has been disabled for this content