Development of the Elk LAPW Code

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
Vext(r)	ρ(r)
Bext(r)	m (r)
Eext(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.

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)].

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 α²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)].

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 .