# Development of the Elk LAPW Code

K. Dewhurst, S. Sharma

Elk is an all-electron full-potential linearised augmented plane-wave (LAPW) code for determining the properties of crystalline solids. In our institute the code is actively developed as a testbed for new methods. It is also extensively used for production, especially for materials which are particularly sensitive to the types of approximation used or for which pseudopotential methods are not appropriate. One aspect which is unique to Elk is that almost all features can be used in combination with each other, resulting in powerful and robust code.

The basic aim of Elk is to solve the Kohn-Sham equations for external fields. Every external field has a conjugate density to which it couples. At the moment, there are four external fields available in Elk:

External field |
Conjugate density |

*V*_{ext}(**r**) |
ρ(**r**) |

**B**_{ext}(**r**) |
**m** (**r**) |

**E**_{ext}(**r**)(=-∇*V*(**r**)) |
ρ(**r**) |

**A** |
**j** |

Although these fields should be related by Maxwell's equations, Elk treats them as independent in order to simulate various conditions outside the crystal. The Kohn-Sham equations are solved in two-step process. In the *first-variational* step a Hamiltonian containing only the scalar potential and **E** field is constructed

and diagonalised with

.

Then in the

*second-variational* step, the magnetic fields, spin-orbit coupling
and

**A** field are added using the first-variational step as a basis

Instead of the usual approach of separating the Kohn-Sham equations into spin-up
and spin-down orbitals, densities and potentials, Elk treats magnetism as

*non-collinear* for which the basic variables are the scalar density

*ρ*(

**r**)
and the magnetisation vector field

**m**(

**r**).
This approach requires fully variational spinor wavefunctions

for which the density and magnetisation are given by

and

respectively.

*Spin-spiral* states also allowed:

in which case the magnetisation density is incommensurate with the crystal
unit cell, having the form
**m**(**r**)=[m_{x}(**r**)sin(**q**·**r**),
m_{y}(**r**)cos(**q**·**r**),m_{z}(**r**)].
This spiral state may be lower in energy than the usual ferromagnetic or
anti-ferromagnetic states for certain materials such as *γ*-Fe
[*Phys. Rev. B* **66**, 14447 (2002)].

#### The LAPW basis

The LAPW basis is constructed by partitioning space into spheres around atoms,
called muffin-tins, and the remaining interstitial region. Linear combinations
of atomic-like orbitals make up the basis functions in the muffin-tins and
plane waves fill the interstitial region. At the
boundary, the functions, and possibly their derivatives, are matched ensuring continuous
or differentiable basis functions. These augmented plane wave functions are
not orthogonal, requiring that the Kohn-Sham equations be solved as a
generalised eigenvalue equation:
*H❘Φ*_{i} ⟩ =* O❘Φ*_{i}〉
where
*O*_{pq} = 〈B_{p}❘B_{q}〉
is the overlap matrix consisting of
inner products of LAPW basis functions |*Bsub>prang;*.
The energy of the radial part of the
muffin-tin functions is chosen to be close to a Kohn-Sham eigenvalue, making
this basis one of the most efficient possible for solids.
- Figure 1: Schematic of the LAPW basis. The basis functions consist of a plane wave augmented by a linear combination of atomic-like functions in the muffin-tins, chosen to so that the resulting function is continuous across the boundary.

An important aspect for any DFT code is the ability to
solve Poisson's equation efficiently for the relevant boundary conditions.
This is needed in order to determine the Hartree
potential but, in the case of Elk, it is also used with complex charge densities
for computing Coulomb matrix elements of the form

It is essential then that the algorithm be as fast as possible, and
considerable effort has been expended in writing a highly
optimised version of Weinert's method [J. Math. Phys.

**22**, 2433 (1981)]
for solving Poisson's equation in the LAPW basis.
These Coulomb matrix elements are used in some of the more advanced features of
the code like Hartree-Fock, the optimised effective potential (OEP) method,
reduced density matrix functional theory (RDMFT),
time-dependent density functional theory (TDDFT), and the Bethe-Salpeter
equation (BSE).

#### Exchange-correlation functionals

Several local density approximation (LDA) and generalised gradient approximation
(GGA) functionals are implimented natively in Elk. It is also
possible to link to the ETSF libxc functional library
(

http://www.tddft.org/programs/octopus/wiki/index.php/Libxc
) which contains
almost every LDA and GGA functional yet conceived. The ability to use
meta-GGA (mGGA) functionals in conjunction with libxc has also been added.
This type of functional requires the kinetic energy density

in addition to

*ρ* and

**m**.
In the case of non-collinear magnetism, Elk employs a technique,
[J. Appl. Phys.

**63**, 3482 (1988)], in which the local 2×2 spin
density matrix

diagonalised and the eigenvalues are used as up and down spin densities in
a LSDA functional. The inverse operation is applied to the
up and down potentials to produce a 2×2 matrix

from which is extracted
the exchange-correlation magnetic field. It is less clear how GGA functionals
should be used for non-collinear situations, but in Elk's case the spin density is
first locally diagonalised into up and down, and then the gradients are computed.
A drawback to these approaches is that the exchange-correlation magnetic field
B

_{xc}(

**r**) is always collinear with the magnetisation density. Thus from the
Landau-Lifshitz-Gilbert equation

we find that for LSDA in the absence of external fields

*∂***m**(

**r**,

*t*)/

*∂**t* = 0,
meaning that there can be no spin dynamics.
This problem can be fixed by use of intrinsically non-collinear functional
approximations which depend explicitly on |

**m**(

**r**)| and its
gradients, or by use of the OEP method for spin-polarised systems as
implimented in Elk [

*Phys. Rev. Lett.* **98**, 196405 (2007)].

#### DFT+U and spin-current tensor moments

Elk has the most sophisticated implimentation of DFT+U available. It can be
used in conjunction with spin-orbit coupling, non-collinear magnetism and
spin-spirals, and has the ability to interpolate between around mean field (AMF)
and the fully localised limit (FLL) [*Phys. Rev. B* **67**, 153106 (2003)].
- Figure 2: Cover of the September 2009 issue of Physical Review Letters showing the spin-current density matrix for two components of a hidden order parameter - the rank five multipole tensor
**w**^{615} - in URu_{2}Si_{2} generated by Elk [F. Cricchio, F. Bultmark, O. Grånäs and L. Nordström, *Phys. Rev. Lett.* **103**, 107202 (2009)].

Recently, the ability to decompose the spin-current density tensor was added
to the code. This allows for the discovery of previously hidden order parameters
in solids which are not simply localised spin moments, but rather specific
moments of the full spin-current tensor.

#### Developments completed at the MPI

- Reduced density matrix functional theory (RDMFT) for solids has been
implimented in Elk. A simple functional was devised and parameterised by
fitting to the exchange-correlation energy of the electron gas. This was used
to determine the band gaps of strongy-correlated insulators.
- Phonon dispersions and electron-phonon coupling matrix elements.
- Calculation of the Eliashberg function α
^{2}*F*(ω)
and solving the Eliashberg equations.
- Inverse random phase approximation (RPA) dielectric function ε
^{-1}
(**G**,**G'**,**q**,w) available as causal, time-ordered
or Matsubara.
- Non-linear optical (NLO) second-harmonic generation.
- Generating and solving the Bethe-Salpeter equation (BSE).
- Linear-response time-dependent density functional theory (TDDFT).

Three forms of parallelism are implemented in Elk, and all can be used in combination with each other, with efficiency depending on the particular task, crystal structure and computer system. OpenMP works for symmetric multiprocessors, i.e. computers that have many cores with the same unified memory accessible to each. The message passing interface (MPI) is particularly suitable for running Elk across multiple nodes of a cluster, with scaling to hundreds of processors possible. Phonon calculations utilize a simple form of parallelism with the code simply examining which dynamical matrices are missing from the run directory: this allows many computers to work independently on a single phonon calculation.

#### Current and future developments at the MPI

- Continuing development of TDDFT including a new functional: the
`boot-strap' kernel, which is able to predict excitonic peaks in optical spectra.
This is also being evaluated for the
*q > 0* case.
- Real-time evolution of solids under the influence of ultra-short laser pulses
in the attosecond or low femtosecond range.
This is a flagship project of our institute and involves new methods for performing
unitary time evolution of the Kohn-Sham orbitals as well as the development of
time-dependent functionals suitable for condensed matter.
- New superconducting density functional theory (SCDFT) method for
*ab-initio* prediction of the superconducting critical temperature *T*_{c}.
- Implementation of magnetic linear response functions for non-collinear
systems in conjunction with DFT+
*U* for direct comparison with experiment as
well as for use in the SCDFT functionals.
- Calculation the spectral density function for solids using RDMFT.
- Development of intrinsically non-collinear functionals beyond the usual
method of diagonalising the spin-density matrix and LSDA.
- An implimentation of the all-electron GW approximation, as well as extending
this to include higher order terms in the Hedin equations
[
*Phys. Rev.* **139**, A796 (1965)].

- Figure 3: Simulated STM image of lithium atoms on a graphene sheet using the local density of states (LDOS) function. Diagonalisation of the first-variational state was performed with a parallel iterative algorithm. The unit cell contains 216 atoms with the lithium atoms highlighted.

We also have begun a biennial workshop, the CECAM Elk Tutorial, aimed at introducing new users to the code and the methods used. The inaugural Tutorial took place in Lausanne from 18-23 July 2011 where experts in various techniques used by Elk were invited to give presentations on their subject. These were recorded and are available online at http://elk.sourceforge.net/cecam.html .

In the spirit of scientific openness, the code is released under the GNU General Public License allowing for the free modification, scrutiny and redistribution of the code. Elk can be downloaded from elk.sourceforge.net .

To top