Spectral quadrature for the first principles study of crystal defects: Application to magnesium

Let $\Omega=L_{1}\times L_{2}\times L_{3}$ be a cuboidal computational domain aligned with the coordinate axis, and let it be discretized using a regular $n_{1}\times n_{2}\times n_{3}$ grid so that the grid spacing is $h_{i}$ in the $i^{th}$ coordinate where $L_{i}=n_{i}h_{i}$. We index the grid points with $q$ and let $K_{\Omega}$ be the set of all grid points. We discretize the gradient and Laplace operators using finite differences on this grid [17, 16].

In each iteration of SCF, we need to determine the spectral nodes and weights at each grid point. As evident from (2.2) and (29), this requires the calculation of the action of Hamiltonian on vectors $v_{k}^{q}$. These vectors are compactly supported around a ball centered at the $q^{th}$ grid point. Therefore, it is sufficient to work with the Hamiltonian $\mathcal{H}^{q}$ projected on to a smaller subspace around the $q^{th}$ grid point. Specifically, we work with a cuboidal region of side $2R_{cut}$ and centered at the $q^{th}$ grid point which we call region of influence. This is shown schematically in Fig. 2a. The controllable parameter $R_{cut}$ depends on the order of the quadrature which in turn depends on material properties and electronic smearing temperature $\theta$ [43, 55].

For grid points close to the edge of $\Omega$, the nodal Hamiltonian can extend beyond the computational domain for grid points which are near the boundary of $\Omega$. However, since we consider problems where our computational domain is embedded in an infinite crystal, we have to compute the Hamiltonian on an extended region $\Omega^{\prime}$ of size $(L_{1}+2R_{cut})\times(L_{2}+2R_{cut})\times(L_{3}+2R_{cut})$, and centered at the origin. The values of the Hamiltonian associated with the annular $\Omega^{\prime}\setminus\Omega$ is obtained using precomputed electronic fields ($\rho$ and $\phi$) for the perfect crystal.

We use domain decomposition for parallel implementation. The computational domain is partitioned into disjoint domains, $\Omega=\bigcup\limits_{p=1}^{N_{p}}\Omega_{p}$, where $\Omega_{p}$ denotes the domain local to the $p^{th}$ processor, and $N_{p}$ is the total number of processors. The collection of grid points belonging to the $p^{th}$ processor is denoted by $K_{\Omega}^{p}$, where $K_{\Omega}=\bigcup\limits_{p=1}^{N_{p}}K_{\Omega}^{p}$, and $K_{\Omega}^{p}\bigcap K_{\Omega}^{q}=\emptyset$ (null set) for $p\neq q$. In our implementation, we use the DMDA class available through the Portable, Extensible Toolkit for Scientific Computation (PETSc) [3, 4] for mesh management. The communication between processes is handled by Message Passing Interface (MPI) libraries [19, 20].

The region of influence of an grid point $q$ owned by a process $p$ (i.e. $q\in K_{\Omega}^{p}$) may extend to the spatial regions owned by neighboring processes. In such a situation, the values of the effective potential $V_{eff}$ from neighboring processes are communicated to the process $p$ using an MPI communicator. In Fig. 2b we schematically illustrate the parallel communication pattern involved for communicating the effective potential $V_{eff}$ from neighboring processes. This reduces the number of MPI related calls otherwise required for matrix vector products.

In each SCF iteration, the spectral weights $\{w_{k}^{q}\}_{k=1}^{K}$ and nodes $\{\lambda_{k}^{q}\}_{k=1}^{K}$ are computed at every grid point $q\in K_{\Omega}$ from the projected Hamiltonian $\mathcal{H}^{q}$. Overall, this is computationally the most expensive part of the method. However, the computation is local and with no MPI related calls. Further, we do not explicitly store the matrices, and their multiplication with a vector is performed in a matrix-free manner. These lead to excellent parallel efficiency.

Once the spectral weights $\{w_{k}^{q}\}_{k=1}^{K}$ and spectral nodes $\{\lambda_{k}^{q}\}_{k=1}^{K}$ are computed at all the grid points $q\in K_{\Omega}^{p}$ for all the processes $p$, we first solve the following equation for the Fermi energy $\lambda_{f}$:

We utilize Brent’s algorithm [7] for this purpose. The summation across the polynomial degree $k$ and the grid point $q$ is performed locally in each processes, and the summation across the processes $p$ is performed with one global MPI communication call.

Once the Fermi energy is calculated, we calculate the point-wise band structure energy $u^{q}$, the point-wise entropy $s^{q}$ and the point-wise electron density $\rho^{q}$, following (17), (18), (19). These are all local.

Once the electron density $\rho$ is calculated at the grid points, we calculate the total electrostatic potential $\phi$ at the grid points by solving the Poisson equation (9) on $\Omega$ subject to Dirichlet boundary conditions obtained from the perfect crystal outside using the generalized minimal residual algorithm (GMRES) [47]. Once the electrostatic potential is calculated at every grid point $q\in\Omega$, we calculate the effective potential $V_{eff}$ at every grid point

where $V_{xc}(\rho^{q})$ is the exchange correlation potential.

The convergence of the SCF iteration is accelerated by employing mixing schemes. In every SCF iteration, we mix the effective potential $V_{eff}$, where we have the option of employing Anderson mixing [2], Pulay mixing scheme [44] and its periodic [5] and restarted [42] variants.

We check the convergence of SCF iteration by calculating the normalized error in the effective potential.

The free energy is computed once the SCF iteration has converged using a discrete version of Eqn. 11

where $\varepsilon_{xc}^{q}$ is the contribution of the exchange correlation energy at the $q^{th}$ grid point, and $E_{self}^{h}$ is the discrete representation of the self energy (3):

In evaluating them, the sum over the grid points are carried out locally on each MPI process followed by a global sum across all the MPI processes $p$.

The final step is the computation of the atomic forces (Eqn. 12). This has two parts, the contributions of the local pseudopotential and the non-local pseudopotentials. The contribution of the local pseudopotential to the atomic force is calculated as

where $\nabla_{h}$ is the gradient operator in the discrete setting. The summation over the the grid points $q$ is local to every process, followed by one summation over the MPI processes.

The non-local contribution to the atomic force is calculated by employing the Clenshaw-Curtis quadrature described in Section 2.2. At each grid point $q\in K_{\Omega}$, we calculate the discrete Chebyshev vectors ${\bf v}_{q}^{k}$ using the iterative scheme (29), and the Clenshaw-Curtis quadrature weights $c_{q}^{k}$ using (28) with the discrete nodal Hamiltonian ${\bf H}_{q}$. The non-local force is given by

The calculation of the atomic forces scales linearly with the number of atoms.

We first verify the convergence with respect to spectral quadrature. We take the same degree $K$ for both the Gauss and the Clenshaw-Curtis quadrature. We study HCP Mg and BCC MgAl. For each of these systems, we randomly perturb the atoms from the ground state to obtain the test configuration, and use a mesh of $h=0.6$ Bohr for pure Mg and $h=0.5$ Bohr for MgAl.

Fig. 3a shows the convergence of the energy and Fig. 3b shows the convergence of atomic forces with $K$. The error is computed with respect to the reference that is taken to be the solutions for $K$=$240$. The decay is not monotone because neither the free energy functional (Eqn. 1) nor the atomic forces (Eqn. 12) is variational with respect to $K$. However, it broadly follows an exponential decay (Error $\approx Ae^{-\beta K}$). The best fit $\beta$ for the various cases is shown in Table 1. The pure Mg system has a higher rate of convergence than the MgAl system. Furthermore, in both the cases, the rate of convergence in total energy and the local component of the atomic force is similar, whereas the rate of convergence in the total force, which require calculating the non-local force component is smaller by almost a factor 2. This is so because the local component of the atomic force and energy use only Gauss quadrature which depends on Lagrange polynomials whereas the non-local component of atomic forces depends on Clenshaw Curtis quadrature with Chebyshev polynomials. The former is known to be more accurate than the latter [52].

Since we require an accuracy of about $1\times 10^{-3}$ and $1\times 10^{-4}$ Ha/atom in our energy, we need between $60$ and $80$ for Mg, and between $60$ and $90$ for MgAl. Similarly, since we require an accuracy between $1\times 10^{-3}$ and $1\times 10^{-4}$ Ha/Bohr in the total force, we need between $50$ and $100$ for Mg and MgAl. These guide our further calculations.

A loose upper bound of the size $R$ of Lanczos vectors in (2.2) is given by $2R\leq hn_{o}K$, where $h$ is the average mesh spacing, and $n_{o}$ is the finite difference order. This suggests a choice of $2R_{cut}=hn_{o}K$. Using the choice of parameters used for this study, $R_{cut}\approx 288$ Bohr. This would make the calculations extremely expensive. However, we now show that it is not necessary by studying the error as a function of $R_{cut}$.

For the Mg and and MgAl systems, we use a Lanczos polynomial degree $K_{L}$ = $80$ and a Chebyshev polynomial degree $K_{C}$ = $100$, and vary $R_{cut}$ from $4.0$ to $14.0$ Bohr. Fig. 4a and 4b shows the error in energy and atomic force as a function of $R_{cut}$. The reference values of energy and atomic force is calculated using an $R_{cut}$ = $16$ Bohr. From these figures, we observe that as $R_{cut}$ increases, the error in energy and forces decrease. We once again observe generally exponential decay, though it is not monotone (also see [43]). This allows us to treat $R_{cut}$ as a controllable approximation parameter.

Next, we turn to the scaling and efficiency of the formulation and parallel implementation developed in this work. We choose magnesium crystal with one aluminum solute atom, placed at the center of the computational domain. The simulation parameters used are $h$ = $0.6$ Bohr, $R_{cut}$ = $12.0$ Bohr, $K_{L}$ = $80$ and $K_{C}$ = $100$ with a desired accuracy of $2\times 10^{-4}$ Ha/atom in energy and $2\times 10^{-4}$ Ha/Bohr in atomic force.

We first perform a strong scaling study with a $600$ atom system, varying the number of cores from $30$ to $420$. The wall times for each SCF iteration is presented in Fig 5a. The parallel efficiency on $420$ relative to 30 cores is $92.6$ percent. This data is plotted in terms of speed up in Fig. 5b. Next, we perform a weak scaling study by selecting systems with sizes ranging from $64$ atoms to $2560$ atoms, while increasing the number of processors from $30$ to $900$. We choose these such that the number of atoms per MPI process is between two and three. Fig 5c shows that the variation of CPU time for one SCF iteration versus the number of atoms is linear ($\approx\mathcal{O}(N)$). Fig. 5d shows that the memory required also scales linearly with respect to the number of atoms. All of this shows excellent scaling of our algorithm. This is due in part to restricting much of the parallel communication to the MPI processes that are neighbors, and keeping the global communication to a minimum.

We note that the spectral quadrature step accounts for greater than $98$ percent of the total time in each SCF iteration. The prefactor of spectral quadrature can be significantly reduced by incorporating reduced basis methods such as Discrete Discontinious Basis Projection ([62]). Further, the number of SCF iterations to achieve a fixed target SCF error increases with system size in metallic systems due to charge sloshing [24]. The introduction of real space preconditioning schemes [50, 29] is likely to reduce this for large metallic systems.

We calculate the formation energy of a monovacancy for various computational domain size from $63$ atoms to $1151$ atoms, and the results are shown in Fig. 6a. We observe that the formation energy strongly depends on cell size, and converges at cell sizes of approximately $1000$ atoms to $0.8456$ eV. This is broadly in agreement with values reported in the literature: calculated values $0.779$ - $0.768$ eV using local pseudo-potential and a coarse grained approach with $1024$ to $1$ billion atoms [40] and measured values of $0.58$ - $0.90$ eV [31, 22, 58, 60]. Fig. 6b shows the electron density on the basal plane in the vicinity of the vacancy. Unsurprisingly, it is depleted at the vacancy. An interesting feature is that the electron density does not display the reflection symmetry of the basal plane (e.g. about the red dashed line). This is because the three dimensional crystal breaks this symmetry at the plane at a height $c/2$ above and below this plane. This is emphasized by indicating the atoms on this upper and lower planes with a ‘*’ in the figure. This observation plays a role in binding.

Next we compute the dilute impurity energy of an aluminum solute atom, and this is shown in Fig. 6c. We again observe that this energy depends on the cell size and converges at a few hundred atoms. The electron density on the basal plane is shown Fig. 6d. The electron density if elevated near the Al solute and the distribution is more symmetric.

We now study the interaction between a pair of point defects: vacancy-vacancy, vacancy-solute and solute-solute. Fig. 7a, shows the six possible nearest neighbor positions on a HCP lattice where one defect is placed at the site marked 0 while the other defects are placed in the sites marked $1$ through $6$ denote the first six nearest neighbor configurations.

Table 2 shows the computed binding energy (cf. (37)) for various defect pairs. Recall that a positive binding energy indicates an energetic propensity for the defects to bind. We observe that all three types of defect pairs – vacancy pair, vacancy-solute and solute-pair – have an energetic preference for binding. This is especially strong for vacancy-solute pairs and solute pairs. Further the binding energy increases with decreasing inter-defect distance in all three cases, indicating that the binding is stronger with decreasing distance. These results are consistent with the observation that voids form readily in Mg, Al has a high diffusivity in Mg and that Al-rich precipitates from readily in Mg. Our results agree quantitatively with the only case (first nearest-neighbor vacancy-solute-pair) where we are aware of experimental observation [45] – Table 2.

The table also shows that while our results agree with those of Ponga, Bhattacharya and Ortiz [40] in the case of vacany pairs, they do not agree with other prior DFT calculations. Specifically, our results predict much stronger binding compared to other results ([59] for vacancy pairs, [1, 51] for vacancy-solute, and [25, 30] for solute pairs). In fact, previous results [25, 30] for solute pairs also predict negative binding energies in some instances.

The previous calculations [59, 1, 51, 25, 30] all used relatively small computatonal cells (64-192 atoms) compared to the current results (1151 atoms). Fig. 7 shows the binding energy of various defect pairs computed for various sizes of the computational domain. We see a significant dependence of the binding energy on the size of the computational cell with convergence only at around a 1000 atoms. Further, there is no common trend. For a solute pair, small cells predict a lack of binding while large cells predict binding. The vacancy-vacancy binding energy is large for small computational domains, but decreases with domain size. For vacancy-solute, we see that the binding energy is small for small computational domains, but increases with domain size. This tells us that the disturbance in the electronic and atomistic fields decay slowly. Further, the previous studies use periodic boundary conditions which means that the defects interact if the unit cell is too small. We have verified that a periodic calculations with a 96 atom unit cell (using SPARC[16] which uses a similar finite difference formulation as the current work) reproduces the result of Uesugi et al. [59] who used the same pseudopotential.

We also note that Ponga, Bhattacharya and Ortiz [40] used a local pseudopotential. In other words, the role of the pseudopotential is less important in these defects compared to that of the slow decay.

To understand the nature of the binding and explore further the decay, we examine the difference in electron density between that of a defect pair with that of a pair of defects:

where $\rho^{0}$ is the electron density in a system with a defect placed at site 0, $\rho^{n}$ is that with a defect at site $n$, $\rho_{p}$ is that of the perfect crystal, and $\rho^{0n}$ is that of a defect pair (we need $\rho_{p}$ to make the system charge neutral). These are plotted for second nearest-neighbor pair (i.e, a basal pair) in Fig. 8. We see that the core of the vacancy pair is severely depleted of electron density while the center of the solute pair has excess electron density. The vacany-solute complex increases electronic density and depletes it at the distal region of the solute. Further, we see that oscillations spread over several atomic spacings before decaying. Furthermore, we see that the electron density in the case of vacancy and solute pairs are not symmetric with respect to the axis through the defects. This is because of the influence of atoms on the basal plane just above and below the plane of observations. Together, these mean that we need large computational cells for accurately computing the binding energies to accurately account for this decay.