Initial demonstration implementation of Hutter Solver

Milestone 2 - This will comprise much of the framework of a releasable implementation of TDDFT using Hutter's published method. Issues of data distribution and parallel efficiency will not be addressed at this stage and an `off-the-shelf' eigensolver will be used. Initial testing and debugging. Code will be benchmarked against a set of existing calculations on molecular systems chosen from scientific literature. This work comprises stages 3-5 of the work plan.

Hutter reformulated the equations of time-dependent Hartree-Fock (TDHF) theory (applied to TDDFT with pure density functionals) such that they could be efficiently implemented in a plane-wave basis set. The TDHF equations are a non-Hermitian eigenvalue equation, which gives the excitation energies directly

$$\left( \begin{array}{ll} \mathbf{A} & \mathbf{B} \\ \mat... ...egin{array}{c} \mathbf{X} \\ \mathbf{Y} \end{array}\right).$$ (3)

In the Tamm-Dancoff approximation, the $\mathbf{B}$ matrices are set to zero, which leads to a Hermitian eigenvalue problem $\mathbf{AX} = \omega\mathbf{X}$. It is this Hermitian eigenvalue problem that will be implemented for the ``Hutter solver''.

In the context of a plane-wave code, the vectors $\mathbf{X}$ and $\mathbf{Y}$ are the plane-wave coefficients of electronic wavefunctions, i.e. $\vert \Phi_{i}^{\left( \pm \right)} \rangle$ from above. In CASTEP, these are four-dimensional arrays, the dimensions being G-vectors (plane-waves), k-points, bands and spin. For a TDDFT calculation, one typically has $\sim 10^5-10^6$ G-vectors, $\sim 100-1000$ bands, a single k-point and spin of 1 or 2 dpending on whether the system is being treated as spin-degenerate. Of course, storing and diagonalising a matrix of this size is innapropriate, especially as one is only interested in the lowest few excited states. Instead we turn to iterative methods, where only the result of the operator on a vector is required.

Hutter defines this operator in equation 35 of his paper, where he chooses to separate it into two components, such that $\mathbf{A} = \mathcal{A} + \mathcal{B}$. The first contribution is from Kohn-Sham orbital energy differences

$$\mathcal{A}\left\vert\Phi^{(-)}_i\right\rangle = (H^{(0)} - \epsilon_i)\left\vert\Phi^{(-)}_i\right\rangle.$$ (4)

The action of this operator was already implemented for DFPT in CASTEP. The second contribution is

$$\mathcal{B}\left\vert\Phi^{(-)}_i\right\rangle = P_c \delta V_{\textsc{scf}}[n^{(-)}]\left\vert\Phi^{(0)}_i\right\rangle,$$ (5)

where $\delta V_{\textsc{scf}}[n^{(-)}]$ is the self-consistent reponse potential for a change in electron density

$$n^{(-)}(\mathbf{r}) = \sum^N_{i=1} \Phi_i^{(0)*}(\mathbf{r})\Phi_i^{(-)}(\mathbf{r}),$$ (6)

and

$$\delta V_{\textsc{scf}}[n^{(-)}] = \int d\mathbf{r^{\prime}... ... n(\mathbf{r^{\prime})}} \right\rbrace n(\mathbf{r^{\prime}}),$$ (7)

which is the same as used in DFPT, but for a factor of 2 in the response density.

Given the effect of the $\mathbf{A}$ operator on a given wavefunction, an `off-the-shelf' iterative eigensolver with a reverse communication interface can be used. By using a library routine, we can thoroughly test our implementation of the operator. We used two different solvers, namely ARPACK³ and EA19, the latter being from the HSL (formerly the Harwell Subroutine Library)⁴. ARPACK is written in Fortran77 and implements the Arnoldi process. HSL-EA19 is written in F95 + TR 15581 and implements a Jacobi-conjugate preconditioned gradients scheme.