{'fraction_non_alphanumeric': 0.02903011934078424, 'fraction_numerical': 0.012762833870051147, 'mean_word_length': 4.634071504802561, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large scale ab-initio simulations of dislocations

We present a novel methodology to compute relaxed dislocations core configurations, and their energies in crystalline metallic materials using large-scale \emph{ab-intio} simulations. The approach is based on MacroDFT, a coarse-grained density functional theory method that accurately computes the electronic structure but with sub-linear scaling resulting in a tremendous reduction in cost. Due to its implementation in \emph{real-space}, MacroDFT has the ability to harness petascale resources to study materials and alloys through accurate \emph{ab-initio} calculations. Thus, the proposed methodology can be used to investigate dislocation cores and other defects where long range elastic defects play an important role, such as in dislocation cores, grain boundaries and near precipitates in crystalline materials. We demonstrate the method by computing the relaxed dislocation cores in prismatic dislocation loops and dislocation segments in magnesium (Mg). We also study the interaction energy with a line of Aluminum (Al) solutes. Our simulations elucidate the essential coupling between the quantum mechanical aspects of the dislocation core and the long range elastic fields that they generate. In particular, our quantum mechanical simulations are able to describe the logarithmic divergence of the energy in the far field as is known from classical elastic theory. In order to reach such scaling, the number of atoms in the simulation cell has to be exceedingly large, and cannot be achieved with the state-of-the-art density functional theory implementations.


Introduction
Dislocations are the main carrier of plasticity in crystalline materials. Their behavior and interactions have strong consequences in the mechanical response of materials. As such, dislocations have been extensively studied since the seminal work of Volterra [1]. Dislocations have long been analyzed by separating their effect into two main contributors, the dislocation core and its longrange field [2]. The former is characterized by strong non-linear quantum mechanical effects; the latter is a long-range elastic field characterized by smooth atomic displacements that are typically modeled using continuum elasticity theory in either its isotropic or anisotropic versions. However, the elastic fields are long-range: the interatomic distances or strains decay inversely with distance ∼ r −1 , where r is the distance from the defect, so that the energy density is logarithmically divergent. This means that the core and far-field are intimately coupled and the traditional separation is problematic. This has many physical implications including the effect of the far-field stresses on the core and the interaction between the dislocation and solute atoms. Therefore, a full understanding requires ab-initio techniques [3]. However, this is not a trivial task since the slow decay of the far field and the intimate coupling means that a very large number of atoms are needed to properly describe this behavior. Unfortunately, this cannot easily be achieved with existing ab-initio methods, since they scale poorly (O(N 3 e ) in traditional approaches and O(N e ) in newer linear scaling approaches [4,5]). Different multiscale approaches have been proposed to overcome these difficulties. Most of these models consist of patching together heterogeneous models at different scales, from DFT to continuum elasticity to couple the long-range effect and the dislocations core. The coupling between the various models varies from parameter passing to hybrid Hamiltonians [6,7,8,9,10,11]. Since these methods bring together distinct models that embody different physics and differing mathematical formulations, they typically assume separation of scales (as in parameter passing) or additional physics or constraints at the interface (as in hybrid Hamiltonians). Often times, the models need to be adjusted or calibrated to the particular phenomenon under consideration, which ultimately detracts from the predictiveness and fundamental nature of first-principles calculations. A variation is the interesting work of Woodward and collaborators [12,13] where they used full density functional theory in a core region and a linear response (Green's functions) outside. Still, the approach separates the domain into multiple regions with two distinct theories.
In this work, we propose a new approach where the equations of density functional theory are solved on a domain large enough to capture both the core and elastic fields in a completely seamless manner. It builds on the recently developed Coarse-Grained Density Functional Theory (CG-DFT) technique developed by Ponga et. al. [14,15] and has it roots in the linear scaling method of Suryanarayana et al. [16]. The latter reformulates the original Kohn-Sham equations [17,18] into density matrix form, and then approximates the density matrix using spectral Gauss quadrature rules. This enables the calculation of the electron density and other quantities of interest at any point in space at fixed or O(1) cost. Consequently the algorithm has linear scaling. However, the fact that the method is local allows the construction of an adaptive numerical approximation that has sub-linear scaling. The key idea is to exploit the decay so that we have full resolution near the core where the details are important and only sample the electronic fields far away where the atomistic displacements decay in a smooth fashion. However, this requires care: even though the atomistic displacements decay smoothly enabling a quasi-continuum approximation [19,20], the electronic fields (like the electron density and Hartree potential) oscillate in an (almost) periodic manner on the scale of the lattice. Our key observation is that the electronic fields of interests can be written as a sum of two parts, a predictor which captures the periodicity in the far-field and a corrector that is complex at the core but decays away from it. We compute the former inexpensively via unit cell calculations at representative points, while we compute the latter on a gradually coarsening set of electronic sampling points by exploiting the local nature of our formulation. Together, this gives sub-linear scaling with the number of electrons allowing the multiscale simulations of hundreds of thousands of atoms with DFT. The technique has been implemented in the MacroDFT code. We will refer to the CG-DFT technique as well as the code as MacroDFT.
The key difference between our prior work [15] is the treatment of the predictor field. In our previous work [15], the predictor was taken to coincide with the periodic solution of the perfect stress-free crystal. While this choice is effective for defects that decay quickly, i.e., r −2 or r −3 , it is not for dislocations where the decay is much slower, i.e., r −1 . The slower decay requires us to update the predictor non-uniformly depending on the far-field deformation, and this adds significant complexity to the method. This is accomplished in this work.
We use MacroDFT to compute the formation energy of prismatic dislocation loops and basal dislocations in magnesium (Mg). Dislocations in Mg have been extensively studied using ab-initio techniques [21,22,23], molecular statics [24] or combination of them [25]. However, often times these simulations are carried out using incompatible boundary conditions such as periodic boundary conditions or free surfaces to patch the interface with continuum approaches, that ultimately affects the energy of the dislocation core, and possible its configuration [3]. These conditions do not reflect the true nature of a dislocation core, and can influence the structure of the dislocation core [26,27,28] as in the case of body-centered materials. Due to these artifacts, no simulation technique using ab-initiotechnique has been able to prove the theoretical logarithmic divergence predicted by the elasticity theory. In this work, we endeavor to address some of the limitations from other ab-initio techniques and study dislocation cores using compatible boundary conditions, and to investigate the scaling of the energy for large simulation cells. Prismatic dislocations loops are made by removing a certain amount of atoms from the bulk material in the basal plane, leading to a Burgers vector (b) in the c−component of the hexagonal closed-packed (hcp) structure of Mg. On the other hand, basal and prismatic screw dislocations are characterized by a much shorter Burgers vector, and their far field is characterized by the elastic decay. Thus, we generate a single dislocation in an infinite crystal, we measure the evolution of the strain energy as a function of the distance from the dislocation core. Our simulations show that the strain energy diverges logarithmically with the ratio r/b, in agreement with the elasticity theory of dislocations [2]. We then compute the interaction energy of the dislocation core with a row of solute atom, and show the transition barriers as dislocation moves through it. The manuscript is organized as follows. We start in Section 2 by introducing the methodology in terms of a density matrix formulation and its spectral representation in Riemann-Stieljes integrals. We then introduce the coarse-grained representation. Next, the computational set-up and verification of the method is discussed in Section 3. The results for prismatic loops and screw dislocations are discussed in Section 4, and the parallel performance of the method is demonstrated in 5. Finally, we summarize the work with main outcomes in Section 6.
where ψ n are the electronic orbitals, ρ(x) = 2 N e/2 n |ψ n (x)| 2 is the charge density and the exchangecorrelation is given by with ε c to be as proposed by Perdew and Wang [30] fitted to accurate Monte Carlo simulations carried out by Ceperley et. al. [31]. The expensive Coulombic double sums are sidestepped by introducing the electrostatic potential φ as the solution to the Poisson equation [32,33]. We minimize this functional over all orbitals subject to orthonormality and this leads to the nonlinear eigenvalue problem with λ n the ordered eigenvalues of the Hamiltonian with φ the Coulomb potential obtained by solving the Poisson's equation denotes the total charge density of the nuclei, with b J (x, R J ) representing the regularized charge density of the J th nucleus, and This nonlinear eigenvalue problem is solved using fixed point iteration in the self-consistent formulation (SCF). This leads to an expensive O(N 3 e ) computational procedure that critically restricts the size of the systems to 100 − 1, 000 electrons in practice. Recent real-space implementations using efficient finite element basis and finite differences can increase this number up to roughly 10,000 electrons [34,35,36].
With the goal of obtaining a linear scaling algorithm, we reformulate the problem above using the the density matrix formulation following Anantharaman and Cancès [37] and Wang et al. [38]. The problem of minimizing (1) is equivalent at zero temperature to the variational problem where the maximum is taken over {γ ∈ χ : |∇|γ|∇ ∈ χ, 0 ≤ γ ≤ 1, Tr γ = N }, χ is the set of all trace class operators on L 2 (R 3 ), and E * xc is the Legendre transform of E xc . One can see this formally by setting γ = Ne n=1 |ψ n ψ n |, using the electrostatic potential to rewrite the Coulomb interactions, using a duality transform for the exchange correlation and exchanging the order of extremization. These steps are explained in detail and justified in [38]. It is possible to show that the solution to variational problem is given by the density matrix and λ f is the Fermi level corresponding to the constraint Tr γ = N , the electrostatic potential φ that satisfies the Poisson's equation (5), the exchange correlation potential u that satisfies (6) and the charge density ρ(x) = γ(x, x).
This density matrix formulation is the basis of linear scaling algorithms (O(N e )). Typically, the density matrix (8) is expanded using polynomials and the expansion truncated assuming bandedness of the Hamiltonian and the decay of off-diagonal components.
We follow a different approach, the linear scaling spectral Gauss quadratures (LSSGQ) following Suryanarayana et al. [16]. The basic idea is to appeal to the spectral representation of the operator γ and to approximate the spectral integrals using Gauss quadratures rules. For any function η ∈ L 2 (R 3 ), where σ is the spectrum of H, and µ ηη is the spectral measure of H contracted with η. We treat the integral as a Riemann-Stieljes integral and approximate it with quadratures rules with quadrature nodes λ η k and quadrature weights w η k to obtain the approximation shown. We use a sixth order finite difference approximation with N f grid points corresponding to a uniform grid spacing h, and choose an orthonormal basis {η p } corresponding to this finite difference basis to represent the Hamiltonian and other quantities. Then, it is possible to show that the ground state energy is given by [15,16] with It then remains to determine the quadrature points λ ηp k and weights w ηp k . We do so using Lanczos iteration [39] (see [16,15] for details). The key observation is that since the finite difference basis has finite support, we can compute the the quantities λ p k , w p k with O(1) effort at each nodal point p. So, the method automatically has linear scaling. Further, this observation means that the evaluation of all the electronic quantities of interest can be done locally (once the Fermi level is known), and this enables the coarse-graining approach described below.
Once we compute the electronic fields, we can compute the forces on the atomic nuclei (the variation of the total energy with respect to the nuclear positions R) easily using the Hellman-Feynman theorem.

Coarse-grained extension
Our goal is a sublinear scaling algorithm that will enable computation of large domains as necessary for dislocations. We observe that the interatomic distances (or strain) decay at the rate of r −1 as we go away from the dislocation core, and the corresponding electronic fields (electron density, electrostatic potential) decay to their periodic behavior at the same rate. So we introduce a  Figure 1: Schematic of an electronic field and its coarse grained representation in MacroDFT. a) Electronic fields around a vacancy defect in a crystalline material. The surface shows the evolution of the field in a 2D plane of the crystalline structure. b) A two-dimensional contour plot of the electronic field (ρ). We can clearly see the effect of the vacancy (illustrated with an 'x') in the field. c) Same electronic field for the pristine crystal without defects (ρ 0 ). d) Illustration of the corrector field (ρ c = ρ − ρ 0 ) for the vacancy example in a two-dimensional slice. Only values above a cutoff has been plotted to ease the view. e) Evolution of the electronic field across the line A-A'. It is evident that the corrector field has large fluctuations near the defect but changes smoothly far away from it. The corrector has been shifted down to make the visualization better. computational basis that makes it possible to describe all details close to the core where it is necessary, but sample it far away where it is not. Importantly, the approach is seamless and adaptive. Given the different behavior, the coarse graining of the atomic positions and electronic fields are done differently.

Coarse grained representation of atoms
The atomistic distances may vary significantly close to the defect, but decay in a polynomial manner. So, we introduce a subset P a of atoms, called representative atoms (RepAtoms). This subset is dense, i.e., it contains every atom close to the defect but gradually coarsen away from it. We track the positions of these atoms and infer the position of other (non-representative atoms) from the position of the representative atoms by introducing linear finite elements shape functions based on a mesh T a with nodes P a as in the atomistic quasi-continuum formulations [15,19,20,40,41,42,43]. Note that this mesh is Lagrangian. Details of the atomic mesh and visual representation of them are provided in Section 3.3.

Coarse grained representation of electronic fields
The electronic fields are complex close to the defect but decay to a periodic field away from it, since the far-field is almost periodic with a period related to the unit cell. So a smooth interpolation is not appropriate. We get around it representing the electronic fields as a sum of a predictor and a corrector. The predictor is a slowly varying almost periodic function away from the defect and captures the periodicity of the electronic field in the far field. The difference between the actual electronic fields and the predictor, therefore smoothly decays away from defect. This is represented by the corrector. Since the corrector decays smoothly, we can use a smooth interpolation away from the defect. This is shown schematically in Figure 1 with a vacancy as the defect. Figure 1(a) shows the evolution of an electronic field around the vacancy centered in the middle of the picture. For reference, a line A−A' is shown to plot the charge density later on. Next, Figures 1(b) and (c) show a contour plot of the charge density for the vacancy and the pristine sample, respectively. We intuitively see that these fields are very similar with the exception of near the defect (marked with an "x"), where large fluctuations are seen. Thus, if we represent the difference between the charge density for the vacancy and the pristine cell, we obtain the contour map shown in Figure  1(d). We see that only near the vacancy this field is large, and very quickly decays to zero. This is the corrector field, and the pristine solution is called in this case the predictor. Putting everything together, we see in Figure 1(e) the three fields, i) the full solution with the vacancy (ρ) in black; ii) the predictor field (ρ 0 ) in red; and iii) the corrector (ρ c ) in blue 2 . Now, it is possible to reconstruct the full solution at any point by using an interpolation scheme provided by the FE mesh. This is schematically shown with the blue markers in the blue line of Figure 1(e), where we see a dense sampling near the vacancy, but a coarsening as one moves away from it. This is the principle for coarse-grained description of MacroDFT.
To be precise, we introduce a fine uniform spatial grid P f and use representation for the electrostatic field and electron densities at each of these points where φ 0 p , ρ 0 p are the predictors and φ c p , ρ c p are the correctors. We defer the description of the predictor till the next section. To obtain the corrector with sublinear scaling, we compute it only a subset of electronic sampling points (ESPs) P c . We then use a finite element mesh T c over the ESPs and use it to interpolate these fields to the fine mesh: The ESPs are dense close to the defect, but gradually coarsen away from it. Note that this mesh is spatial or Eulerian. Finally, note that we have two meshes, one atomistic and Lagrangian that is used to describe the relaxation of the atoms as the forces and energy are minimized; and an electronic mesh that Eulerian to solve the Kohn-Sham equations in the sample. Both meshes have adaptive resolution, as previously described. They interact through the theory and therefore, it is important to introduce two special functions that map the ESP into the elements of the T QC , and positions to the elements of the fine mesh T f . For more details, the reader is referred to [15].
Total energy 2 The corrector has been shifted down to make the visualization better.

Algorithm 1 MacroDFT algorithm
initialize. Initial guess for the positions of RepAtoms P a and electronic fields at ESPs P c .
repeat Use the atomistic mesh T a to find position of all atoms, Generate the predictor fields φ 0 p , ρ 0 p for each p ∈ P f , repeat Find the corrector fields φ c p , ρ c p using (14) for each p ∈ P f , Find {λ q k } and {w q k } using Lanczos algorithm for each q ∈ P c , Find N e , ρ q for each p ∈ P f using (11), (12), until convergence of the electronic fields, Compute the forces on the RepAtoms P a , Update the position of the RepAtoms P a , until convergence of the atomic positions.
We compute the atomic positions and electronic fields according to the Algorithm 1. It remains to calculate the total energy (10) of the system. We do so by sampling using cluster as follows. Let C q be the set or cluster of nodal points on the fine mesh contained in a ball of radius r around the q th ESP. We approximate the total energy (10) as [15] E ≈ q∈Pc n q E q (15) where n q is the weight of the q th ESP and it is defined as n q = 1 if q belongs to the full resolution zone #P f associated with qth ESP otherwise (16) and F q is the average Free-Energy at the q th ESP computed as

Selection of the predictor field
The key for an efficient and accurate coarse-grained implementation is the choice of the predictors. In our previous work [15], the predictor field was taken to coincide with the periodic solution of the perfect stress-free crystal. While this choice is effective for defects that decay quickly, i.e., r −2 or r −3 , where r is the distance from the defect, it is not for dislocations where the decay is much slower, i.e., r −1 . Therefore, we need a more accurate approximation that reflects the deformation.
Recall the triangulation T QC . In each element of this triangulation, the positions of the atoms is obtained from a linear interpolation and thus periodic with a unit cell deformed by the deformation gradient given by the atomistic element. It is important that the predictor reflect this periodicity. Therefore, we extract the deformation gradient from each atomistic element, deform the unit cell with this deformation gradient, and compute the electronic field associated with this distorted unit cell. We then extend this periodic field to the entire element. Note that this field automatically reflects the local strain. However, it can be discontinuous at the element boundaries; we remedy this by smoothing with a L 2 → H 1 projection. Thus, once the atomic positions are computed with the non-linear conjugate gradient, a new deformation gradient, F e , is computed per each element. The deformation gradient is then applied to a unit cell, and the electronic fields are computed on it. We then extract the corrector fields, i.e., ρ 0 (x) and φ 0 (x), and map them back to the full simulation as schematically shown in Figure  2. The calculations in the unit cell are carried out by individual nodes, and their results are stored on the disk to save memory with an MPI-file. While this puts some restrictions to the speed of the code, the time is not comparable with the amount taken to solve the DFT equations, and then it is affordable.

Computational set-up
We now describe the details of the simulations. The details of the meshes used in this work are described in Table 2.

Verification of the method
Verification of our method has been provided in our previous publication [15], and here we reproduce some of these results for completeness. In our formulation, we first start by investigating two main parameters that condition the accuracy and convergence of the methodology. These two parameters are the number of quadrature points, used to compute Eqs. 10 and 11, and the number of ESP used to discretize the atomic unit cell in full resolution. Figure 3 shows the convergence with respect to number of quadrature points and spatial discretization for a unit cell of Mg. We see how the difference in the energy is reduced as both quadrature points and number of ESPs is increased. In this work, we retained the spatial resolution we used in our previous work. This is about 2,300 ESPs per unit cell, which ensures a convergence of about 0.01 meV per unit cell. We have also used a sufficiently large number of quadrature points, K = 40. With these parameters, is possible to achieve cohesive energy for Mg that are converged to less than a 1 meV. Using a unit cell and (a) Convergence with respect to quadrature points (b) Convergence with respect to spatial discretization  These values are shown in Table 1 for completeness. We can see that the values obtained with MacroDFT are in close agreement with ABINIT calculations. The stable stacking fault energy for Mg was also computed with MacroDFT using the described parameter in a simulation cell containing 2 × 2 × 48 unit cells, i.e., 192 atoms. The values of the stable stacking fault energy was found to be ∼ 27 mJ·m −2 , in good agreement with previous works [25]. Convergence of our coarse-grained simulations has been done in our previous work [15]. In addition, MacroDFT has been used to compute the energy of twin boundaries in Mg, showing good agreement with other techniques [44]. Here, we used the same level of coarsening, that is, the element size increases in different regions of coarse representation by a factor of two. This ensures that the evaluated energy does not depend of the mesh size.

Prismatic dislocation loops in magnesium
A single crystal Mg was generated using the minimum energy configuration for the bulk crystal at a = 3.109Å and c/a = 1.626 [15]. The overall length of the computational domain was 40a 0 × 40a 0 × 40c 0 containing N = 256, 000 atoms. A full resolution zone of 8a 0 × 8a 0 × 8c 0 was provided in the center of the simulation cell. This area provided fully resolved atomic and electronic mesh containing around ∼ 1, 126, 400 ESPs. Surrounding the full resolution zone, multiple coarse-grained regions were provided for both the electronic and atomic fields. Two Delaunay triangulations over the ESPs and RepAtoms were performed, and a FE mesh was provided over these triangulations to interpolate atomic positions and corrector fields. From the center of the full resolution area, N v atoms were removed to generate a cluster of vacancies. The boundary conditions were set up such that the solution decay to the pristine solution at the end of the simulation cell, for both atomic positions and electronic fields. Let u(R) be the atomic displacement field, then u(R) = 0, ∀ R / ∈ P QC and ρ c (x) = 0, φ c (x) = 0 ∀ x / ∈ P f . This means that we enforced the displacement of the atoms to be zero out of the simulation box and forced the corrector field of the electron density and electrostatic potential to be zero outside the simulation cell too. The selection of these boundary conditions is justified based on the idea that the deviations of the electronic field and atomic displacements decay to zero far away of the defect. [1120]{1010} screw dislocations were generated by displacing the atoms in the simulation cell using the isotropic elastic solution [2]. The dislocations were generated such that the dislocation line was aligned with the x− direction in the sample as shown schematically in Figure 4. A full resolution zone of a 0 × 16 √ 3a 0 × 4c 0 surrounding the dislocation core was provided for the a 0 3

Screw dislocations in magnesium
[1120]{0001}. By contrast, a full resolution zone of a 0 × 4 √ 3a 0 × 16c 0 surrounding the dislocation core was provided for the a 0 3 [1120]{1010}. In the full resolution region, the spatial resolution for the electronic degree of freedom is about 2, 300 ESPs per unit cell. The total number of ESPs in the full resolution area is around 281, 600. Surrounding this area, three coarse-grained regions were provided. Each coarse-grained region had a constant element size that increased with the distance to the core. The grand total number of electronic degrees of freedom was around 560,395 ESPs. Over this set of points, a Delaunay triangulation was performed and a FE mesh was generated by using four-nodes tetrahedral elements with linear interpolation functions.
Similarly to the electronic mesh, an atomic mesh was also provided. The atomic mesh also contained a full resolution area that coincides with the full resolution electronic region. Away from the dislocation core, multiple coarse-grained regions were provided. The perfect atomic positions are used as starting point for this mesh. Figure 5(a) shows spatial distribution of the atoms for a a 0 3 [1120]{1010} screw dislocation as well as the FE mesh for a sample with three coarse-grained regions. Figure 5  The coarse grained atoms on the slip plane need to have a special capability -similar to the discontinuous FE Galerkin formulation-where the elements are allowed to have a dual position in order to interpolate properly the atomic positions in both sides of the slip plane. With the screw dislocation generated, a FE mesh was constructed over the representative atoms using an α-Delaunay triangulation. The α-Delaunay triangulation was needed because the set of representative atoms with the screw dislocation is not a convex hull and highly distorted elements were generated with a regular Delaunay triangulation. A full resolution simulation of the previously described dislocation will involve 304 million ESPs. This leads to a total acceleration due to the spatial coarse-grained approach of 545 times.
For all simulations we used a k B θ = 0.8 which is consistent with other DFT studies of metals and the energy and forces were minimized by using the Polak-Ribière version of the non-linear conjugate gradient (NLCG). The boundary conditions were set up such that the atomic displacements outside the computational cell were given by the elastic solution [2]. By doing so, we can compute for all elements outside the computational box the deformation gradient F e a priori. Thus, for each atom outside the computational domain, the electronic field and atomic displacements are fully determined. This means that ρ c (x) = 0, φ c (x) = 0 ∀ x / ∈ P f .

Prismatic dislocation loops in magnesium
The energy of the computational cell is minimized using a NLCG and the formation energy and the binding energy were computed as and where E f Nv is the formation energy for the cluster with N v vacancies, E(N − N v ) is the total energy of the system containing N v vacancies, Nv is the biding energy for N v vacancies, and E f 1 is the formation energy of a single vacancy. Table 3 shows the values of the formation energy and the binding energy for both the noncollapsed and collapsed configurations. The binding energy for cluster of vacancies lying on the basal plane is quite large and correspond to a very stable defect. The binding energy computed using Eq. 19 measures the difference between the energy in a cluster of vacancies and the energy of the same number of individual vacancies. Therefore, such a large binding energy indicates that vacancies will always try to regroup in clusters instead of being randomly disperse in the sample, at least at low temperatures where entropic effects are not important. Figure 6 shows the electron density through two planes at y = 0 and z = 0, respectively for a sample containing 7 vacancies in the non-collapsed configuration.
It is also interesting to analyze the binding energies of the same clusters collapsed, in order to understand the formation of prismatic dislocation loops in Mg. For instance, for the sample containing 7 vacancies, the collapsed configuration was always unstable and the atoms came back to the non-collapsed configuration shown in Figure 6. On the other hand, when the number of vacancies was increased to N v = 19 and N v = 37, the binding energy experienced a change, as shown in Table 3. The binding energy for N v = 19 was around 2.99 eV, and for N v = 27 was around 10.13 eV. This highly non-linear demeanor of the binding energy indicates a preference towards large prismatic dislocation loops and therefore, sets a minimum size of PDL that can be observed in Mg. Although configurations with larger number of vacancies are not feasible with MacroDFT, there is a clear trend in the binding energy which favors large PDLs, as observed in experiments [45,46,47] of Mg and nanoporous Mg alloys.

Screw dislocations in magnesium
After relaxing the dislocation core, the excess of energy for the dislocations was measured with respect to the pristine Mg configuration. To measure the energy of the dislocation core, we took a cylinder of radius R centered in the dislocation core (y = 0, z = 0) with principal axis along the dislocation line. Then, we analyzed the excess of energy for different values of R. The excess energy was defined as where N is the total number of atoms in the cylinder of radius R, E D (N, R) is the total energy in the same cylinder after the core has been relaxed, and E C is the cohesive energy of pristine Mg. The energy in the cylinder, E D (N, R) is computed by evaluating the electronic fields for each ESP inside the cylinder using the predictor-corrector approach. Once the fields are evaluated, the total energy was computed using Eq. 17.
In the analysis of dislocations it is usually common to show the evolution of the excess energy in the cylinder as a function of the ratio R r 0 , where r 0 is an arbitrary constant. It is common practice to take r 0 = b, and under this assumption the energy contained in the cylinder can be expressed as [2] where µ = 19.9 GPa is the shear modulus of Mg computed with MacroDFT, and E core is an indeterminate constant that needs to be evaluated using discrete techniques such as ab-initio or  molecular dynamics calculations. Remarkably, the value E core is not a physical quantity as it depends on the value r 0 that is arbitrary. The logarithmic divergence of Eq. 21 indicates that in order to understand the energy of the dislocation core, the energy has to be sampled across very large samples spanning multiple Burgers vectors in order to obtain the elastic pre-factor µb 2 4π . Figure 9 shows the evolution of the excess energy, computed with Eq. 20, as a function of the log(R/r 0 ). By adjusting the E core values to the results obtained with MacroDFT, we have concluded that E core = 7 meV/Å and E core = 7.5 meV/Å for the basal and prismatic slip systems, respectively. This small strain energy per unit length is very difficult to measure and illustrate the importance of coarse-grained techniques such as MacroDFT. The results also show an excellent agreement for the far field with the elastic solution (characterized by the slope of the blue curve) which lends credence to the prediction of dislocation far field obtained with MacroDFT and therefore, to the coarse-grained scheme used. Figure 10 shows the differential displacement plots [48] where SF stands for stacking fault. By examining the atomic positions as shown in Figure 10(a) we observed that the dissociation distance was around ∼ 2.9a 0 . By comparing with previous studies of screw dislocations in the basal plane, we found that this value is slightly larger than previously reported [25]. This could be attributed to a number of reasons, including different boundary conditions, pseudo-potential used, etc. The a 0 3 [1120]{1010} screw dislocation is stable and does not split into two partials. This is consistent with stacking fault calculations along the prismatic plane. These calculations show that there is not stable local minima in this direction, and therefore, the prismatic dislocation cannot split.  The atoms are shown in green, while the arrows indicate the magnitude of the differential displacement field. The prismatic dislocation has a very compact core, while the core in the basal plane is slighter larger due to the existence of stable stacking fault in Mg.

Aluminum solute-dislocation interaction
We now proceed to measure the interaction energy between an Al solute atom and the dislocations for various position along the slip plane. In order to do so, we replace a Mg atom by an Al solute stating from the position y = 0, z = 0. We only replace atoms lying on the slip plane for the two screw dislocations considered in this work. Then, the excess energy was computed as where E D (N, R, r Al ) is the energy in a cylinder of radius R = 500Å, containing N-1 Mg atoms and one Al atom at r Al , and E D (N, R, r Al max ) is the energy in a cylinder of radius R = 500Å, containing N-1 Mg atoms and one Al atom at the r Al max , which is the maximum distance from the core at which we have placed the solute atom. Figure 11 shows the excess of energy for various positions of the solute atom for the basal and prismatic slip system (Figure 11(a) and (b), respectively). Figure 11 depicts a non-linear interaction energy between solute and dislocation core, and this interaction is very strong near the core decaying very fast relatively quickly. Remarkably, the maximum interaction energy is of the order of 50-60 meV/Angs. which is larger than the core energy of the dislocations.

Parallel performance
We now proceed to discuss the parallel performance and speed-up achieved with MacroDFT in MIRA an IBM BG/Q 1.6 GHz PowerPC A2 supercomputer property of Argonne National Laboratory. Measuring the efficiency and performance are essential for new DFT methodologies and codes, as most of the classical plane-wave codes scale poorly for more than 64 processors. Thus, the next generation of computational tools should be able to harness petascale resources that can further accelerate the discovery of new materials and alloys through accurate ab-initio calculations. This is challenging since both the methodology and implementation need to be highly efficient to handle thousands of processors at the same time. Here, we provide results for the   Each node has 16 physical processors and capable of running 64 threads per node. strong scaling of MacroDFT using up to 262,144 threads. The IBM BG/Q 1.6 GHz PowerPC A2 node architecture has 16 physical processors, each node is capable of running 4 threads. Thus, the user can specify any combination of threads and processor per node such that the number # processors × # threads = 64. The ratio processors/threads has to be evaluated for each code/problem. In our case, we found that the best combination was to use processes with 64 threads, which reduces the internode communication. Using this combination, the threads share the same memory, the number of MPI communications was reduced significantly. For the benchmark, we have used the non-collapsed cluster of vacancies containing 7 vacancies described in Section 3.2 and compute one SCF evaluation using different numbers of threads. Figure 12 shows the parallel performance for MacroDFT from 8,192 (2,048) to 262,144 (65,536) threads (processors). The results show that MacroDFT has a good efficiency up to 131,072 threads with an efficiency of ∼ 87.5%, while for 262,144 threads, an efficiency of ∼ 68.75% is achieved. This remarkable performance is due to the fact that the DFT formulation (Eq. 9) can be solved locally and the relevant electronic fields ESP can be evaluated with only the information of a small neighborhood surrounding the point. Since the computational representation of the ESPs is in realspace, the cost of computing the electronic fields can be distributed between multiple processors taking advantage of massively parallel computers, as only the local components of the Hamiltonian need to be updated every SCF iteration. According to the scaling performance shown in Figure 12, the time taken by the interpolation scheme is negligible in comparison with the time consumed to determine the nodes and weights of the spectral quadrature rule. By contrast, when the number of threads is increased significantly, 262,144 or more, the performance goes down due to multiple factors including load imbalance, few number of ESP per node, interpolation, and data writing to disk. Nonetheless, the performance of MacroDFT is outstanding when compared with the stateof-the-art DFT codes, which cannot even run in such a large number of processors. Therefore, real-space DFT implementations have an opportunity to take advantage of supercomputing power and over-perform traditional plane-wave based implementations.

Conclusions
We have developed a novel approach for the study of dislocations in crystalline materials using ab-initio techniques. In this approach, the equations of density functional theory are solved on a domain large enough to capture both the core and elastic fields in a completely seamless manner. It technique is based on three-main pillars i) reformulation of the Kohn-Sham equations into a density matrix formulation and subsequent spectral representation of the density matrix (See Section 2.1); ii) Separation of the electronic fields in predictor and corrector counterparts (Equations (13), (14) and Section 2.3); and iii) Interpolation scheme for electronic fields and displacements using finite element meshes (2.2). The proposed approach has been implemented in the MacroDFT code [42] which can be executed in MPI, thread and hybrid options according to the architecture of the high performance cluster available. The parallel efficiency of the code has been studied, achieving good parallel performance in extreme scale computers capable of running ∼ 260, 000 threads.
The code and method has been used to study the relaxed configurations of the atoms in multiple dislocation cores in Mg, a crystalline metallic material. For PDLs, our simulations have predicted that vacancies always bind together, at least for low temperatures where entropic effects are weak. In screw dislocations gliding in the basal and prismatic planes, we have found that the strain energy stored in the computational sample scales accordingly to the elasticity theory for both cases. The results indicate that in order to obtain good prefactors, the simulations has to be sufficiently large, of the order of a few hundred thousand atoms. The agreement of the elastic prefactor between our simulations and theory lends credence to the approach, and provides numeric evidence that the coarse-grained approach can be used in complex scenarios where quantum mechanical interactions are important, such as in dislocation cores in alloys, grain and twin boundaries.