Real-time adaptive finite element solution of time-dependent Kohn–Sham equation

Abstract

In our previous paper (Bao et al., 2012 [1]), a general framework of using adaptive finite element methods to solve the Kohn–Sham equation has been presented. This work is concerned with solving the time-dependent Kohn–Sham equations. The numerical methods are studied in the time domain, which can be employed to explain both the linear and the nonlinear effects. A Crank–Nicolson scheme and linear finite element space are employed for the temporal and spatial discretizations, respectively. To resolve the trouble regions in the time-dependent simulations, a heuristic error indicator is introduced for the mesh adaptive methods. An algebraic multigrid solver is developed to efficiently solve the complex-valued system derived from the semi-implicit scheme. A mask function is employed to remove or reduce the boundary reflection of the wavefunction. The effectiveness of our method is verified by numerical simulations for both linear and nonlinear phenomena, in which the effectiveness of the mesh adaptive methods is clearly demonstrated.

Introduction

In 1984, Runge and Gross [24] discovered the one-to-one correspondence between time-dependent one-body density $n(\overrightarrow{x},t)$ and time-dependent one-body potentials $v_{\mathit{ext}}(\overrightarrow{x},t)$ in a system under certain quite general conditions, for a given initial state. Based on Runge–Gross theorem, the time-dependent electron density becomes a fundamental variable of the system, which can generate all other physical properties of the system. Analogous to the Kohn–Sham system, the time-dependent electron density can be obtained by solving a fictitious non-interacting system, the time-dependent Kohn–Sham system, instead of the original interacting one. Consequently, the dimension of the problem may be dramatically reduced, and the theoretical and numerical study of large electronic system becomes possible.

Since then, a great deal of studies have been conducted on developing numerical methods for TDKS equation. In the perturbative regime, the time-dependent linear response theory (TDLRT) is one of the most popular methods to study the response of the system to an external potential $\delta v_{\mathit{ext}}({\overrightarrow{x}}^{\prime },t^{\prime })$. For the method, once the linear density response function $\chi (\overrightarrow{x},{\overrightarrow{x}}^{\prime },t-t^{\prime })$ is obtained, the first-order response of all properties derivable from the density with respect to any scalar field can be calculated smoothly. Such linear density response function can be derived straightforwardly from the ground-state of the system. Hence, the TDLRT is widely used in the calculations of photoabsorption spectra, the excited state, etc., of the given electronic structures. In the cases that the higher-order components in the external potential can no longer be omitted, the time-dependent higher-order response function can be introduced with an analogous derivation to the linear case. Similar to the TDLRT, this higher-order density response function gives the desired high-order response of all properties derivable from the density. To date, the TDLRT and related higher-order methods haven been widely used in the computational chemistry community, and have been realized in a lot of software such as Siesta [28], GPAW [21], [34], and Octopus [5]. However, the drawback of the TDLRT is also quite apparent. From the derivation of the linear and higher-order density response functions, we know that the ground-state perturbation of the system is required. This means that the variation of the time-dependent electron density cannot be very far away from the ground-state electron density, or else the reliability would be destroyed. This drawback limits the applications of TDLRT on the study of the non-perturbative regime of a given system where the time-dependent electron density could vary dramatically in a large domain. Another issue about the TDLRT, from the computational point of view, is that a high number of unoccupied orbitals needs to be calculated to guarantee the accuracy and reliability of TDLRT, which increases the difficulty of the convergence on the ground-state calculation, especially for a large and complex electronic structure system.

Developing numerical methods for real-time time-dependent Kohn–Sham equation (RT-TDKS) can resolve the above mentioned issues of TDLRT effectively. In a RT-TDKS simulation, the system is evolved in the time domain by solving TDKS equation, and the time-dependent electron density $n(\overrightarrow{x},t)$ is recorded at each time instant. Then all time-dependent properties of the system can be derived from $n(\overrightarrow{x},t)$, according to the Runge–Gross theorem. From the description above, we can see that there is no restriction on the time-dependent external potential $v_{\mathit{ext}}(\overrightarrow{x},t)$, and the item that we need to record, say, the time-dependent electron density $n(\overrightarrow{x},t)$, only depends on the occupied orbitals. Furthermore, both linear effect such as photoabsorption spectra and nonlinear effect such as high harmonic generation can be numerically studied by using the above process. Therefore, it appears that RT-TDKS methods are more powerful than TDLRT methods. However, there are several challenges on developing efficient numerical methods for RT-TDKS equation. From the model point of view, so far, an effective and reliable time-dependent exchange-correlation potential $v_{xc}[n](\overrightarrow{x},t)$ is still unavailable. Theoretically, this quantity is unknown analytically, and is global in both time and space. In practice, the approximations such as the adiabatic local density approximation (ALDA), are used in the simulation. Although ALDA is local in both time and space, and has good performance in the perturbative simulations, its performance in the non-perturbative simulations is unacceptable, especially when the external field is sufficiently strong. Hence, to develop a generally effective local approximation of $v_{\mathit{ext}}[n](\overrightarrow{x},t)$ is crucial for the study of RT-TDKS equation. On the other hand, from the computational point of view, more computational resources are needed by RT-TDKS methods compared with TDLRT methods. Unlike TDLRT methods which only focus on an interested frequency of the response of the system, RT-TDKS methods can produce the response of the system at all frequencies in a single simulation. It makes RT-TDKS methods more demanding on the memory and CPU time. In addition, since the nonlinear property of the equation, and the fact that the time propagator of TDKS equation is actually a time-ordered exponential operator, to develop efficient numerical methods for RT-TDKS equation is quite nontrivial. In spite of these challenges, the advantages of RT-TDKS methods have already motivating a lot of research on related fields.

There have been a lot of pioneer works on developing numerical methods for RT-TDKS equation. The related algorithms have also been realized in several libraries, such as QBox [14], [26] where the plane-wave method is used, Siesta [28] where linear combination of atomic orbitals (LCAO) basis sets are adopted, GPAW [21], [34] where the grid-based projector-augmented wave method is used, and Octopus [5] where the finite difference method is used. In these libraries, the frameworks of the algorithm for RT-TDKS are well-established, and excellent results are presented from simulating photoabsorption spectra of the electronic structures to calculating nonlinear phenomena such as high harmonic generation and multiphoton ionization. Although these existing methods and libraries have been widely applied, there is still one common potential issue on their practical applications, i.e., the methods are all rely on the regularity of the domain and/or the structured property of the mesh. In addition, various boundary conditions could appear in the simulation with complex configuration. To handle different boundary condition is still a nontrivial challenge for plane-wave method. As one kind of real-space methods, the finite element methods have the potential to deliver more robust numerical solvers on handling various boundary conditions and domain with complex geometry, even compared with finite difference methods mentioned above which is another kind of real-space methods. In addition, since the finite element methods have been widely used in other research fields such as computational fluid dynamics, computational electromagnetics, etc., there have been a lot of mature techniques for developing efficient methods such as a posteriori error estimations, multigrid methods, etc., which can benefit the development of the efficient methods for RT-TDKS equation. People may refer to [23] for the finite element methods in the electronic structure calculations. Although the finite element methods have been developed for the ground-state calculations of the electronic system in [29], [35], [2], [1], [3], [10], [17], etc., little is known about the application of finite element methods for the RT-TDKS. In [27], a finite element method solver for time-dependent and stationary Schrödinger equations is proposed. However, only one-dimensional problem is considered there. In [6], the finite-difference/finite-element time–space discretization is utilized for the TDKS. Although the applications for the three-dimensional carbon nanotube is presented there, only one-dimensional problem is actually solved with the mode approach. The full three-dimensional discretization for the TDKS equation may be very complicated and CPU-time demanding in the implementation. However, it is meaningful for the practical applications, especially for the system with a general structure.

For three dimensional simulations, Lehtovaara et al. extend their all-electron finite element solver for the Kohn–Sham equation [17] to the RT-TDKS case [18]. Excellent numerical results can be observed there. To generate an appropriate mesh, an elegant method is introduced in [17], [18], say, the mesh is created by merging highly nonuniform but symmetric atomic meshes to a global, nonuniform, and unstructured mesh. Similar mesh generation methods can also be observed from [29], [25]. The meshes generated by these methods can handle ground-state calculations, and the perturbative simulations very well. However, there is a difficulty for using such fixed mesh in non-perturbative simulations since the time-dependent electron density could change dramatically during the simulation. As a result, a sufficiently dense mesh in a sufficiently large domain around the nuclei might be necessary for a reliable simulation. In some case, such domain could be very large, and it is a dynamic process for the variation of the electron density during the simulation. Hence, from the computational point of view, mesh adaptive methods, including the local refinement and coarsening, offer a better choice rather than a global dense mesh. In the simulation, the mesh density is adjusted dynamically according to the numerical results with the mesh adaptive methods, and computational resource can be potentially saved. Although there have been a lot of works on developing mesh adaptive methods for ground-state calculations [10], [35], [8], rare of them is related to the RT-TDKS simulations.

In this paper, we develop an adaptive finite element method for RT-TDKS simulations, based on our previous work in [1]. There are two main components in a numerical discretization for the TDKS equation, the spatial discretization and the temporal discretization. For the spatial discretization of the TDKS equation, we follow [1] to employ the linear finite element discretization. For the temporal discretization of the TDKS equation, a second-order semi-implicit Crank–Nicolson (SOICN) scheme is adopted, which conserves the unitary and the time-reversal symmetry of the time propagator well. Since the TDKS equation is a complex equation, the linear system derived is actually a complex-valued one. We will study this complex-valued system by separating the wavefunction into the real part and the imaginary part, and the final matrix derived by this way is a large block one. Since this block matrix is nonsymmetric, we use the bi-conjugate gradient (BiCG) method to solve the linear system. Other than the nonlocal part of the pseudopotential in Kleinman–Bylander form, all four submatrices of the large block matrix actually have the same sparsity pattern. Based on this observation, an algebraic multigrid (AMG) method is developed in this paper for efficiently solving the block matrix. The numerical results show that the derived AMG solver works very well. As discussed above, the mesh adaptive method could potentially improve the efficiency of the RT-TDKS simulations. In this paper, we follow our previous work in [1] to use the hierarchical geometry tree (HGT) for the mesh data, and introduce a heuristic a posteriori error estimator to handle the mesh refinement and coarsening. With the mesh adaptive methods, meshes with a large diameter can be used close to the boundary of the domain. Consequently, a sufficiently large computational domain can be used for a RT-TDKS simulation, without significantly increasing the number of grid points. Furthermore, our solution will not be very sensitive to the usage of the absorbing boundary conditions. In our simulation, the simple mask function method [19] serves as the absorbing boundary condition.

This paper is arranged as follows. In the next section, the TDKS model is briefly summarized. Then the details of the numerical discretization of the TDKS equation are presented in Section 3. We discuss several numerical issues in Section 4. In Section 5, a variety of numerical simulations are presented to show the effectiveness and reliability of our numerical method. Finally, we conclude the research and discuss some future directions.