Recent Advances in the Korringa-Kohn-Rostoker Green Function Method

The Korringa-Kohn-Rostoker (KKR) Green function (GF) method is a technique for all-electron full-potential density-functional calculations. Similar to the historical Wigner-Seitz cellular method, the KKR-GF method uses a partitioning of space into atomic Wigner-Seitz cells. However, the numerically demanding wave-function matching at the cell boundaries is avoided by use of an integral equation formalism based on the concept of reference Green functions. The advantage of this formalism will be illustrated by the recent progress made for very large systems with thousands of inequivalent atoms and for very accurate calculations of atomic forces and total energies.


Introduction
Since the early days of quantum mechanics it has been recognized that the quantum mechanical laws could be used to describe "a large part of physics and the whole of chemistry" as, for instance, Paul Dirac [1] argues in his paper on "Quantum Mechanics of Many-Electron Systems" in 1929.Therein he also argues that this leads to "equations much too complicated to be soluble" and that "approximate practical methods ... should be developed" which can be used "without too much computation".Within Hohenberg-Kohn-Sham density-functional theory the approximation consists in the expressions chosen for the unknown exchange-correlation functional and the simplification that arises is the use of the electronic density as the fundamental quantity instead of many-electron wave function.But even with this enormous simplification, accurate density computations for solids are not trivial so that many different techniques have been suggested and developed over the decades.
Most techniques presently in use are based on expanding the one-electron Kohn-Sham wavefunction, from which the density is calculated, into a set of judiciously chosen basis-functions and use the minimum principle for the total energy to determine the expansion coefficients.A quite different approach was pursued by Eugene Wigner and Frederick Seitz [2] in their work "On the Constitution of Metallic Sodium" published in 1933.Here, the one-electron Schrödinger equation is solved by dividing the periodic crystal into non-overlapping cells centred at each atom, by determining single-cell solutions and putting them together with appropriate boundary conditions.It is the aim of the present contribution to relate this historical Wigner-Seitz approach to the modern KKR-GF method and to illustrate the capabilities of this approach by showing that very large systems with thousands atoms can be studied and that very accurate forces and total energies can be calculated.

Historical Perspective
The essential point in the Wigner-Seitz cellular method is the idea that the solution in each cell can be obtained individually under the constraint of correctly specified boundary conditions at the surface of the cell.The boundary conditions are particularly simple in periodic crystals if the cells are constructed by planes which bisect perpendicularly the lines connecting the atomic positions.For this construction, which leads to the so-called Wigner-Seitz cells, the wave-function derivative in the direction perpendicular to the surface must vanish because of symmetry.In Ref. [2], Wigner and Seitz further simplify the computations by approximating the cells by "Wigner-Seitz" spheres and by choosing a spherical potential in each sphere.
For many years the Wigner-Seitz method was considered as "perhaps the most highly developed and widely applied technique for obtaining one-electron wave functions in solids" as John Reitz [3] states in his "Solid State Physics" review article "Methods of the One-Electron Theory of Solids" in 1955.About two decades later in 1971, however, John Ziman [4] argued in his "Solid State Physics" review article "The Calculation of Bloch Functions" that the symmetry argument which is "at the heart of the cellular method" and "played in important part in the early history of band-structure calculations" has "outlived its practical usefulness".A plausible reason for this view was probably the observation that the assumption of "sphericity" of the cell and the potential is only meaningful for close-packed solids and that this assumption can be exploited in the KKR method [5,6] with much less computational effort. 1 In the KKR method the numerically demanding wave-function matching at the cell surfaces is avoided by using multiple-scattering theory which goes back to Lord Rayleigh's article [9] "On the Influence of Obstacles arranged in Rectangular Order upon the Properties of a Medium" from 1892.The essential point is that the outgoing wave scattered at each cell can be analytically transformed into incoming waves at the other cells if the waves can propagate freely between different scattering events.The free propagation is guaranteed in the historical KKR method by confining the atomic potentials to non-overlapping muffin-tin spheres.Because the multiple-scattering picture cannot be generalized to space-filling potentials and because the numerically effort in the historical KKR method increases for N atoms as N 4 compared to N 3 in basis-set methods, the KKR method was not considered to be competitive for density-functional calculations.
In principle, this unfavourable situation has changed with the work of John Beeby [10] who pointed out "that it is not necessary to know Bloch functions in the lattice to find the density of the electrons" because the density is also easily obtained from the Green function of the one-electron Schrödinger equation.The essential point in Beeby's formulation is the determination of the Green function by an integral equation where the integral kernel consists of a product of the potential and a free-space Green function.The advantage is that the integral equation can be solved analytically and that only a single numerical approximation is required which consists in choosing a maximal angular momentum l max for the expansions of the free-space reference Green function in spherical harmonics Y lm (as explained, for instance, in Ref. [11]).In this way the KKR-GF method uses an efficient and mathematically clear implementation of the cell-based approach that Wigner and Seitz had in mind.

Recent Advances
The solution of the integral equation for the Green function can be written as 1 The fact that the cellular method can be used for the correct form of the cell [7] and even for non-spherical potentials [8] has practically received no attention.In N 2 scaling mode (black curve with circles) the quadratically increasing total computing time is well distributed over the increasing number of processors leading to a wall clock time which increases proportional to N. A similar parallel efficiency is obtained in N scaling mode (blue curve with squares) where the wall clock time only slightly increases with increasing system size.Here, the matrix elements in Eq. ( 1) were set to zero and thus not calculated if n was outside a region consisting of 959 atoms around each atom n.

EPJ Web of Conferences
Here r and r are cell-centred coordinates in the cells around the atoms n and n , r > and r < are the vectors r and r with larger or smaller length, is a parameter over which must be integrated to obtain the density, R n lm and S n lm are regular and irregular scattering solutions for cell n and G nn L L are matrix elements which are related to the ones of free space by a linear matrix equation.

Large systems
If standard methods, for instance, plane-wave expansions based on pseudopotentials are used, densityfunctional calculations for systems with thousands of atoms are very expensive for two reasons.The computing time usually increases as N 3 and the efficient use of supercomputers with hundred of thousands of processors is complicated because a considerable amount of data must be exchanged between processors.Here, the KKR-GF method is more advantageous if a reference system consisting of artificial repulsive potentials is used instead of free space.Then, as described in Ref. [12], the Green function (1) can be calculated independently for each cell n and each value of l and m with high parallel efficiency (see Fig. 1).The computing time increases as N 2 , if no compromise on the accuracy is made, and only as N, if small errors in total energies are tolerated.In this way, the KKR-GF methods was applied to date for large systems with 4096 cells to study the role of vacancies in metal-insulator transitions of the phase-change material GeSb 2 Te 4 [13], with 1024 cells to study the role of vacancy clustering on the superparamagnetism in Gd-doped GaN [14] and with 2048 cells to study the force field pattern around antisite defects in a shape memory alloy [15].

Forces
In the KKR-GF method the force F n on atom n is calculated by where Z n is the nuclear charge and n n core (r) the density of core electrons of atom n.V M (r) is the Madelung potential (defined as the Coulomb potential calculated without the diverging self-interaction arising from the charge Z n ) and V(r) the effective Kohn-Sham potential.For force calculation it is important that accurate potentials are determined near the nuclear sites.This requires to calculate accurate electron densities near these sites which is complicated because of the diverging behaviour of the irregular solutions S n lm (r) for |r| → 0.Recently it was shown [15] how this problem can be overcome by an integral equation formalism based on a subinterval technique for the radial integration.The essential point is that the individual subinterval solutions are combined by analytical expressions so that no numerical wave-function matching is required.3) that are used for the calculation of all quantities but the single-particle energies.For the calculation of the single-particle energies, potential matrix elements up to l, l = 8 are used.The black curve with circles is for Al, the red curve with squares for Cu and the blue curve with diamonds for Pd.Plotted are the deviations from the results obtained with l max = 8.

Total energies
Although the KKR-GF method yields accurate results for many physical properties, the convergence of calculated total energies with respect to l max is often considered to be less satisfactory.The reason for this is that usually not enough potential matrix elements are taken into account because of computational resource limits.Recently, this problem was analysed [11] with the result that for millielectron-volt accuracy (see Fig. 2) only the single-particle part of the total energy must be calculated with high values of l and l .All other quantities in the KKR-GF method can be calculated with low values of l max which substantially reduces the computing time because the single-particle energies must be calculated only once after the self-consistency steps.
a e-mail: ru.zeller@fz-juelich.deDOI: 10.1051/ C Owned by the authors, published by EDP Sciences, This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Figure 1 .
Figure 1.Wall clock time for one density-functional self-consistency step for systems from 1024 to 16384 atoms.In N 2 scaling mode (black curve with circles) the quadratically increasing total computing time is well distributed over the increasing number of processors leading to a wall clock time which increases proportional to N. A similar parallel efficiency is obtained in N scaling mode (blue curve with squares) where the wall clock time only slightly increases with increasing system size.Here, the matrix elements in Eq. (1) were set to zero and thus not calculated if n was outside a region consisting of 959 atoms around each atom n.

Figure 2 .
Figure2.Total energy results as function of the highest angular momentum l max of the potential matrix elements (3) that are used for the calculation of all quantities but the single-particle energies.For the calculation of the single-particle energies, potential matrix elements up to l, l = 8 are used.The black curve with circles is for Al, the red curve with squares for Cu and the blue curve with diamonds for Pd.Plotted are the deviations from the results obtained with l max = 8.