---
abstract: 'A general real-space multigrid algorithm MIKA (Multigrid Instead of the K-spAce) for the self-consistent solution of the Kohn-Sham equations appearing in the state-of-the-art electronic-structure calculations is described. The most important part of the method is the multigrid solver for the Schrödinger equation. Our choice is the Rayleigh quotient multigrid method (RQMG), which applies directly to the minimization of the Rayleigh quotient on the finest level. Very coarse correction grids can be used, because there is in principle no need to represent the states on the coarse levels. The RQMG method is generalized for the simultaneous solution of all the states of the system using a penalty functional to keep the states orthogonal. Special care has been taken to optimize the iterations towards the self-consistency and to run the code in parallel computer architectures. The scheme has been implemented in multiple geometries. We show examples from electronic structure calculations employing nonlocal pseudopotentials and/or the jellium model. The RQMG solver is also applied for the calculation of positron states in solids.'
address: 'Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100, FIN-02015 HUT, FINLAND'
author:
- 'T. Torsti, M. Heiskanen, M. J. Puska, and R. M. Nieminen'
title: 'MIKA: a multigrid-based program package for electronic structure calculations'
---

=10000

[2]{}

Introduction {#sec:introduction}
============

The goal of computational materials science and also that of modeling of nanoscale man-made structures is to calculate from first principles the various chemical and/or physical properties. This requires the solution of the electronic (and ionic) structures of the system in question. The density-functional theory (DFT) [@kohn98] makes a huge step towards this goal by casting the untractable problem of many interacting electrons to that of noninteracting particles under the influence of an effective potential. However, in order to apply DFT in practice one has to resort to approximations for electron exchange and correlation such as the local-density approximation (LDA) or the generalized-gradient approximation (GGA). Moreover, in the case of systems consisting of hundreds or more atoms it is still a challenge to solve numerically efficiently for the ensuing Kohn-Sham equations.

We have developed a real-space multigrid method called MIKA (Multigrid Instead of the K-spAce) for the numerical solution of the Kohn-Sham equations [@mgarticle1]. In real-space methods[@beckrev; @arias; @waghmare], the values of the wave-functions and potentials are presented using three-dimensional point grids, and the partial differential equations are discretized using finite differences. Multigrid methods[@brandt1; @beckrev] overcome the critical slowing-down (CSD) phenomenon occuring with basic real-space relaxation methods. Several approaches employing the multigrid idea have appeared during recent years[@briggs; @ancilotto; @fattebert2; @wang1].

From the different multigrid methods available for the solution of the Schrödinger equation, we have picked up the Rayleigh Quotient Multigrid (RQMG) method introduced by Mandel and McCormick [@McCormick]. This approach differs from full-approximation-storage[@brandt2; @beck1; @wang1; @costiner] (FAS) methods, as well as from those methods[@briggs], where the eigenproblem is linearized.

In the RQMG method the coarse grid relaxation passes are performed so that the Rayleigh quotient calculated on the [*fine*]{} grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid and the coarse grid representation problem is avoided. Mandel and McCormick[@McCormick] introduced the method for the solution of the eigenpair corresponding to the lowest eigenvalue. We have generalized it to the simultaneous solution of a desired number of lowest eigenenergy states by developing a scheme which keeps the eigenstates separated by the use of a penalty functional[@mgarticle1].

Numerical Methods {#sec:methods}
=================

In our RQMG application the coarse grid relaxations are performed by the so-called coordinate relaxation method. It solves the discretized eigenproblem $$H u = \lambda B u$$ by minimizing the Rayleigh quotient $$\label{Ray}
  \frac{\langle u\arrowvert H\arrowvert u\rangle}
       {\langle u\arrowvert B\arrowvert u\rangle}.$$ Above, $H$ and $B$ are matrix operators chosen so that the Schrödinger equation discretized on a real-space point grid with spacing $h$ is satisfied to a chosen order $O(h^n)$. In Eq. (\[Ray\]) $u$ is a vector containing the wave function values at the grid points. In the relaxation method, the current estimate $u$ is replaced by $ u' = u + \alpha d$, where the search vector $d$ is simply chosen to be unity in one grid point and to vanish in all other points, and $\alpha$ is chosen to minimize the Rayleigh quotient. This leads to a simple [^1] quadratic equation for $\alpha$. A complete coordinate relaxation pass is then obtained by performing the minimization at each point in turn and these passes can be repeated until the lowest state is found with desired accuracy.

Naturally, also the coordinate relaxation suffers from CSD because of the use of local information only in updating $u$ in a certain point. In order to avoid it one applies the multigrid idea. In the multigrid scheme by Mandel and McCormick[@McCormick] the crucial point is that [*coarse*]{} grid coordinate relaxation passes are performed so that the Rayleigh quotient calculated on the [*fine*]{} grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid. In practice, a coarse grid search substitutes the fine grid solution by $$\label{rqmgchgeq}
 u_f' = u_f + \alpha I_c^f e_c,$$ where the subscripts $f$ and $c$ stand for the fine and coarse grids, respectively, and $I_c^f$ a prolongation operator interpolating the coarse grid vector to the fine grid. The Rayleigh quotient to be minimized is then $$\begin{aligned}
\label{rqmgeq}
&    \frac{\langle u_f + \alpha I_c^f d_c \arrowvert H_f \arrowvert
         u_f + \alpha I_c^f d_c \rangle}
       {\langle u_f + \alpha I_c^f d_c \arrowvert B_f \arrowvert
         u_f + \alpha I_c^f d_c \rangle} = \qquad \qquad \qquad \qquad \nonumber \\
&\qquad \qquad \qquad   
      \frac{ \langle u_f \arrowvert H_f u_f \rangle
         + 2\alpha \langle I_f^c H_f u_f \arrowvert d_c \rangle
         + \alpha^2 \langle d_c \arrowvert H_c d_c \rangle
       }
       { \langle u_f \arrowvert B_f u_f \rangle
         + 2\alpha \langle I_f^c B_f u_f \arrowvert d_c \rangle
         + \alpha^2 \langle d_c \arrowvert B_c d_c \rangle
       }.\end{aligned}$$ The second form is obtained by relating the coarse grid operators, $H_c$ and $B_c$, with the fine grid ones, $H_f$ and $B_f$, by the Galerkin condition $$\label{galerkincond}
   H_c  = I_f^c H_f I_c^f; \quad
   B_c  = I_f^c B_f I_c^f; \quad
   I_f^c = \left(I_c^f\right)^T.$$ The key point to note is that when $H_f u_f$ and $B_f u_f$ are provided from the fine grid to the coarse grid, the remaining integrals can be calculated on the coarse grid itself. Thus one really applies coordinate relaxation on the coarse grids to minimize the *fine level* Rayleigh quotient. This is a major departure from the earlier methods, which to some extent rely on the ability to represent the solution of some coarse grid equation on the coarse grid itself. Here, on the other hand, one can calculate the *exact* change in the Rayleigh quotient due to *any* coarse grid change, no matter how coarse the grid itself is. There is no equation whose solution would have to be representable.

In the MIKA package we have generalized the RQMG method to the simultaneous solution of several mutually orthogonal eigenpairs. The separation of the different states is divided into two or three subtasks. First, in order to make the coarse grid relaxations converge towards the desired state we apply a penalty functional scheme. Given the current approximations for the $k$ lowest eigenfunctions, the next lowest, $(k+1)$’th state is updated by minimizing the functional $$\label{rqmgneq}
 \frac{\langle u_{k+1}\arrowvert H\arrowvert u_{k+1}\rangle}
      {\langle u_{k+1}\arrowvert B\arrowvert u_{k+1}\rangle}
 + \sum\limits_{i=1}^{k}
    q_i \frac{\left|\langle u_i | u_{k+1}\rangle\right|^2}
             {\langle u_i | u_i\rangle \cdot
              \langle u_{k+1} | u_{k+1}\rangle}.$$ The minimization of this functional is equivalent to imposing the orthonormality constraints against the lower $k$ states, when $q_i \rightarrow \infty$. By increasing the shifts $q_i$ any desired accuracy can be obtained, but in order to obtain a computationally efficient algorithm a reasonable finite value should be used, for example $$q_i = (\lambda_{k+1}-\lambda_i) + {\rm Q},$$ where $Q$ is a sufficiently large positive constant. In our test calculations $Q$ is of the order of $Q=0.5\ldots 2$ Ha.

The substitution (\[rqmgchgeq\]) is introduced in the functional  (\[rqmgneq\]) and the minimization with respect to $\alpha$ leads again to a quadratic equation. This time the coefficients contain terms due to the penalty part.

While the penalty functional keeps the states separated on the coarse levels, we apply a simple relaxation method (Gauss-Seidel) on the finest level. The Gauss-Seidel method converges to the nearest eigenvalue, so ideally no additional orthogonalizations would be needed. In practice, however, we use Gramm-Schmidt orthogonalizations and subspace rotations[@mgarticle1]. However, the number of fine grid orthogonalizations remains quite plausible, for example, in comparison with the conjugate gradient search of eigenpairs employing only the finest grid [@seitsonen].

The Kohn-Sham equations have to be solved self-consistently, [*i.e.*]{} the wave functions solved from the single-particle equation determine via the density (solution of the Poisson equation and the calculation of the exchange-correlation potential) the effective potential for which they should again be solved. To approach this self-consistency requires an optimized strategy so that numerical accuracy of the wave functions and the potential increase in balance, enabling the most efficient convergence [@wang1; @waghmare]. Our strategy in MIKA for self-consistency iterations is illustrated in Fig. \[fig:strategy\]. The Poisson equation for the Coulomb potential is solved also by the multigrid method.

Examples {#sec:results}
========

We have demonstrated[@mgarticle1] the performance of the MIKA scheme in calculating the electronic structures of small molecules and solid-state systems described by pseudopotentials. As a typical application Fig. \[fig:deep\_state\] shows the electron density of the so-called deep state localized at a neutral, ideal vacancy in bulk Si. It was shown, that the accuracy of 1 meV for the total energy was reached after three or four V-cycles, and that the amount of cpu-time needed was of the same order as when applying state-of-the-art plane-wave codes. We obtained an average convergence rate of approximately one decade per self-consistency iteration. This is of the same order as those reported by Wang and Beck [@wang1] in their FAS scheme or by Kresse and Furthmüller [@kresse2] in their plane-wave scheme employing self-consistency iterations. The convergence rate of one decade per self-consistency iteration is better than that obtained by Ancilotto [*et al.*]{} [@ancilotto] in the FMG scheme and much better than the rate reached in the linearized multigrid scheme by Briggs [*et al.*]{} [@briggs].

We have applied the RQMG method also for the calculation of positron states in solids. Fig. \[fig:positrons\] shows how the delocalized positron state in the perfect $\alpha$-quartz is trapped in to a Si-vacancy. Positron states are a particularly simple case for our method, because only the lowest energy wave function needs to be calculated in a given potential, so that no orthogonalizations or penalty functionals are needed. Moreover, in a simple scheme an electron density calculated without the influence of the positron can be used as the starting point [@pos_rev]. However, even for the positron states the superior performance of the multigrid method in comparison with straightforward relaxation schemes is evident. For example, we have calculated the positron state in the Si-vacancy in bulk Si using a supercell containing 1727 atoms. The solution of the positron wave-function using the RQMG-method took less than a minute of cpu-time on a typical work station. To put this in proper context, J. E. Pask [*et al*]{} [@Pask] report a similar calculation, based on the finite element method but without multigrid acceleration, for a supercell containing 4096 Cu atoms. The result converged within 1 ps took ’ just 14.3 hr ’ of CPU time.

We have also applied the MIKA scheme in two-dimensional problems for quantum dots employing the current-spin-density functional theory (CSDFT), see Ref. [@Henri]. Moreover, we have implemented the RQMG-method in cylindrical coordinates enabling very efficient and accurate calculations for atomic chains, or systems which can be described using axisymmetric jellium models. Fig. \[fig:cylinder\] shows a selected wavefunction of a system where a chain of four carbon atoms is sandwhiched between two planar jellium leads.

Summary and outlook {#sec:conclusions}
===================

In the MIKA program package the RQMG method introduced by Mandel and McCormick [@McCormick] is generalized for the simultaneous solution of a desired number of lowest eigenenergy states. The approach can be viewed to belong to a third group of multigrid methods, in addition to FAS and techniques where the eigenproblem is linearized. In principle, one can use arbitrarily coarse grids in RQMG, whereas in the other multigrid methods one has to be able to represent all the states also on the coarsest grid.

We are convinced that our method will compete with the standard plane-wave methods for electronic structure calculations. However, some straightforward programming is still required. Implementation of the Hellmann-Feynman forces, required for the optimization of the ionic structures is under way.

During the RQMG V-cycle, the states are all relaxed simultaneously and independently of each other. A parallelization over states would therefore be natural to implement on a shared memory architecture. We have parallelized the MIKA codes over k-points, and over real-space domains. The domain decomposition is the appropriate method for distributed memory parallel computers.

We acknowledge the contributions by Henri Saarikoski, Paula Havu, Esa Räsänen, Tero Hakala, and Sampsa Riikonen in sharing their experience of the use of the MIKA package in different applications and preparing the figures \[fig:positrons\] (T.H.) and \[fig:cylinder\] (P.H.). T.T. acknowledges financial support by the Vilho, Yrjö and Kalle Väisälä foxundation. This research has been supported by the Academy of Finland through its Centre of Excellence Programme (2000 - 2005).

W. Kohn, Rev. Mod. Phys. [**71**]{}, 1253 (1998). M. Heiskanen, T. Torsti, M. J. Puska, and R. M. Nieminen, Phys. Rev. B [**63**]{}, 245106 (2001). T. L. Beck, Rev. Mod. Phys. [**72**]{}, 1041 (2000). T. A. Arias, Rev. Mod. Phys. [**71**]{}, 267 (1999). U. V. Waghmare, H. Kim, I. J. Park, N. Modine, P. Maragakis, and E. Kaxiras, Comp. Phys. Comm. 137, 341 (2001). A. Brandt, Math. Comp. [**31**]{}, 333 (1977). E. L. Briggs and D. J. Sullivan and J. Bernholc, Phys. Rev. B, [**52**]{}, 5471 (1995); [*ibid.*]{}, [**54**]{}, 14362 (1996). F. Ancilotto, P. Blandin, and F. Toigo, Phys. Rev. B, [**59**]{}, 7868 (1999). J.-L. Fattebert, J. Comput. Phys. [**149**]{}, 75 (1999). J. Wang and T. L. Beck J. Chem. Phys. [**112**]{}, 9223 (2000). J. Mandel and S. McCormick, J. Comput. Phys. [**80**]{} (1989) 442. A. Brandt, S. F. McCormick and J. W. Ruge, SIAM J. Sci. Stat. Comput. [**4**]{}, 244 (1983). T. L. Beck and K. A. Iyer and M. P. Merrick, Int. J. Quantum Chem. [**61**]{}, 341 (1997). S. Costiner and S. Ta’asan, Phys. Rev. E, [**52**]{}, 1181 (1995). A. P. Seitsonen and M. J. Puska and R. M. Nieminen, Phys. Rev. B, [**51**]{}, 14057 (1995). G. Kresse and J. Furthmüller, Comput. Mat. Sci. [**6**]{}, 15 (1996). M. J. Puska and R. M. Nieminen, Rev. Mod. Phys. [**66**]{}, 841 (1994). J.E. Pask, B.M. Klein, P.A. Sterne and C.Y. Fong, Comp. Phys. Comm. [**135**]{}, 1, (2001). H. Saarikoski, M. J. Puska, and R. M. Nieminen, in this conference proceedings.

[^1]: For the sake of simplicity, the wave-functions, and thus $\alpha$, are here assumed real. We have implemented the complex case as well.
