Massively parallel density functional calculations for thousands of atoms: KKRnano

Applications of existing precise electronic-structure methods based on density functional theory are typically limited to the treatment of about 1000 inequivalent atoms, which leaves unresolved many open questions in material science, e.g., on complex defects, interfaces, dislocations, and nanostructures. KKRnano is a new massively parallel linear scaling all-electron density functional algorithm in the framework of the Korringa-Kohn-Rostoker (KKR) Green's-function method. We conceptualized, developed, and optimized KKRnano for large-scale applications of many thousands of atoms without compromising on the precision of a full-potential all-electron method, i.e., it is a method without any shape approximation of the charge density or potential. A key element of the new method is the iterative solution of the sparse linear Dyson equation, which we parallelized atom by atom, across energy points in the complex plane and for each spin degree of freedom using the message passing interface standard, followed by a lower-level OpenMP parallelization. This hybrid four-level parallelization allows for an efficient use of up to $100000$ processors on the latest generation of supercomputers. The iterative solution of the Dyson equation is significantly accelerated, employing preconditioning techniques making use of coarse-graining principles expressed in a block-circulant preconditioner. In this paper, we will describe the important elements of this new algorithm, focusing on the parallelization and preconditioning and showing scaling results for NiPd alloys up to 8192 atoms and $65536$ processors. At the end, we present an order-$N$ algorithm for large-scale simulations of metallic systems, making use of the nearsighted principle of the KKR Green's-function approach by introducing a truncation of the electron scattering to a local cluster of atoms, the size of which is determined by the requested accuracy. By exploiting this algorithm, we show linear scaling calculations of more than $16000$ NiPd atoms.

Figure 1

Reduction of the number of required iterations by application of the third initial guess strategy as described above for the system Ge${}_1$Sb${}_2$Te${}_4$. The reduction ratio $\gamma$ is defined as the ratio of the sum of all iterations needed for all sites, $\mathbf{k}$ and energy points, as well as $lm$ components with and without application of the initial guess. Here, we considered system sizes from $N=64$ to 512 and display $\gamma$ as a function of the degree of convergency of the self-consistency steps, which is represented by the rms error, the variance of actual ($s$), and previous ($s-1$) potential $\propto \parallel V_{\left(s\right)}\left(r\right)-V_{(s-1)}\left(r\right)\parallel$, from which $X_{\left(s\right)}^{\left(0\right)}=X_{(s-1)}$ is taken. The average reduction is plotted as a straight line.

Figure 2

In all five panels, a two-dimensional supercell (boundary marked with gray dashed lines) is drawn and the five different conceptual steps needed for the block-circulant preconditioning scheme are illustrated. (a) Shows the supercell on the level of atoms, which is here exemplified by a two-dimensional disordered lattice of two arbitrary types of atoms (light and dark blue) with overall $N=64$ atoms. In addition, the supercell is partitioned into $M_{\mathrm{bl}}=16$ subblocks, where each subblock contains $N_{\mathrm{bl}}=4$. The borders of the subblock are indicated by blue lines. Then, overall $M_{\mathrm{bl}}^x=4$ and $M_{\mathrm{bl}}^y=4$ subblocks in the $x$ and $y$ directions are required to represent the full supercell. (b) Illustrates how those subblocks are labeled over the entire supercell by the index $j$ from $j=1$ to $j=M_{\mathrm{bl}}=M_{\mathrm{bl}}^x{M}_{\mathrm{bl}}^y$. In addition, it is shown how this index $j$ is related to the position ($x_j$, $y_j$) in real space. In (c), the interblock interactions of subblock $j=6$, which is highlighted by thick blue lines in (b) and (c), are depicted in the space of the row index $i_l\left(j\right)$. For the sake of simplicity, index $i_l\left(j\right)$ is running exclusively over nearest-neighboring subblocks. Further, the relative geometrical position of the interacting blocks ${\Delta }_{i_l\left(j\right)}^x$ and ${\Delta }_{i_l\left(j\right)}^y$ are specified for this simplified example. (d) Schematically shows the full interaction matrix of the supercell from all $M_{\mathrm{bl}}=16$ blocks amongst each other, highlighting two selected types of interaction $i_l\left(j\right)=1$ (orange colors) and $i_l\left(j\right)=2$ (blue colors), where the variations in color represent variations in the individual interactions. In direct contrast, (e) depicts the consequence of averaging the full interaction matrix to effective interactions by means of Eq. (29), which are accordingly represented by uniform colors. Note that (d) and (e) are schematical illustrations being not one-to-one related to (a)–(c).

Figure 3

(a) Convergency of the norm of the residual vector $\left|r\right|$ as a function of the number of TFQMR iterations combining the initial guess approach and block-circulant preconditioning (BCP) for an arbitrary column of a matrix corresponding to Si${}_{432}$. Four different qualities of starting vectors rated by the initial residual $|r^{\left(0\right)}|=A{X}^{\left(0\right)}-B$ are shown for comparison with (thick lines) and without (thin lines) applied preconditioning. Note that a log-log scale is used. In (b) and (c), the CPU time required to obtain convergency is shown for all qualities of initial guesses used in (a) in consistent color coding. Timings with and without applied preconditioning are shown in (c) and (b), respectively. Note that the scale is different in (b) and (c). For the sake of comparison, the CPU time required to set up the BCP is specified in (c). The CPU time to apply the initial guess is for all scenarios negligibly small and not shown.

Figure 4

Sum of TFQMR and TFQMR+BCP iterations over all $lm$ components and one arbitrary atom in the test system Ni${}_5$Pd${}_{251}$. Here, all 27 energy integration points are considered and no initial guess has been used. A random displacement from the ideal atomic sites on the order of 1$\%$ has been introduced on all sites. The inset shows the distribution of the energy points along the integration contour. The color-coded points illustrate the corresponding positions of energy points in the main graph.

Figure 5

Schematic workflow and nested parallelization of four levels of hierarchy as implemented in KKRnano. Starting with a serial process (gray), the computational work is distributed to the real-space parallelization (black). Here, three branches correspond to three atoms. These processes are further split, e.g., into two processes distributing operations over energy integration points (blue). Subsequently, those processes can be further split in spin-up and spin-down processes (red) in case of spin polarization. While up to this point all parallelization has been implemented in a distributed memory approach using MPI 2.0, the inner nested parallelization (orange) of the solution of the Dyson equation is based on OpenMP parallelization, which is in this example split into two threads.

Figure 6

Scaling of KKRnano on a Blue Gene/P (Ref. 25) for a Ni${}_x$Pd${}_{1-x}$ system as a function of number of atoms per unit cell $N$ from $N=$108 to 6912 with a Ni concentration of $x\approx 3\%$ using one $\mathbf{k}$ point and always one processor per atom. Block-circulant preconditioning has been applied with blocks of four atoms, and no initial guess has been used. Displayed are the execution time $t_f$ (blue), the total CPU time $t_{f}N=t_f{N}_{\mathrm{p}}$ (orange), as well as total CPU time without the time required for MPI communication (gray crosses). Straight lines are linear $cN$ (blue) and quadratic $c{N}^2$ (orange) fits to the data points, where $c$ is a constant, which is slightly different for both fits. All time measurements have been executed using the performance analysis tool scalasca (Ref. 27).

Figure 7

Speedup of KKRnano as a function of number of atoms $N$ in the supercell. The number of processors is always equal to the number of atoms $N$ defining the speedup of 1. All optional levels of parallelization are probed separately on a Blue Gene/P (Ref. 25) architecture for Ni${}_x$Pd${}_{1-x}$ alloys. The details are the same as described in the caption of Fig. 6. Labels $\{{\mathrm{p}}_{\mathrm{S}},{\mathrm{p}}_{\mathrm{E}},{\mathrm{t}}_{\mathrm{OMP}}\}$ with ${\mathrm{p}}_{\mathrm{S}}$, ${\mathrm{p}}_{\mathrm{E}}$, and ${\mathrm{t}}_{\mathrm{OMP}}$ specify the number of processes/threads used per spin, energy, and OpenMP parallelization, respectively. For the energy parallelization, dynamic load balancing as explained in the text has been adopted.

Figure 8

Speedup of KKRnano (blue) with respect to the number of processors combining subsequently all levels of parallelization on a Blue Gene/P (Ref. 25) architecture for a Ni${}_x$Pd${}_{1-x}$ alloy of 4096 atoms versus ideal speedup (gray line). The computational details are the same as described in the caption of Fig. 6. Labels $\{{\mathrm{p}}_{\mathrm{S}},{\mathrm{p}}_{\mathrm{E}},{\mathrm{t}}_{\mathrm{OMP}}\}$ with ${\mathrm{p}}_{\mathrm{S}}$, ${\mathrm{p}}_{\mathrm{E}}$, and ${\mathrm{t}}_{\mathrm{OMP}}$ specify the number of processes/threads used per spin, energy, and OpenMP parallelization. For the energy parallelization, dynamic load balancing as explained in the text has been adopted.

Figure 9

Scaling of KKRnano in double-logarithmic representation without (open blue circles) and with applied truncation (filled blue circles). The test system is Ni${}_{1-x}$Pd${}_x$ with $x=5$$\%$ and one $\mathbf{k}$ point on JUGENE. Here, parallelization over atoms, spin, and two groups of energy points are used. The number of processing nodes (4 CPUs per node) is shown in dark gray.