Atom-projected and angular momentum resolved density of states in the ONETEP code

Local and angular momentum projected densities of states (DOS) are invaluable sources of information that can be obtained from density functional theory calculations. In this work, we describe a theoretical framework within ONETEP’s linear-scaling DFT formalism that allows the calculation of local (atom-projected) and angular momentum projected density of states l-p-DOS. We describe four different bases that can be used for projecting the DOS with angular momentum resolution and perform a set of tests to compare them. We validate the results obtained with ONETEP’s l-p-DOS against the plane-wave DFT code CASTEP. Comparable results between ONETEP’s and CASTEP’s charge spilling parameters are observed when we use pseudo-atomic orbitals as the projection basis sets. In general, the charge spilling parameters show remarkably low values for projections using non-contracted spherical waves as the angular momentum resolved basis. We also calculate the d-band and d-band centres for Pt atoms in (1 1 1) facets of cuboctahedral Pt nanoparticles of increasing size, which is an example of l-p-DOS application commonly used as an electronic descriptor in heterogeneous catalysis. Interestingly, the different projection bases lead to similar conclusions, showing the reliability of the implemented method for such studies. The implementation of these methods in a linear-scaling framework such as ONETEP provides another tool for analysing the electronic structure of complex nanostructured materials.

In particular, the local (atom projected) density of states (l-DOS) and the angular momentum projected density of states (p-DOS) are useful tools to obtain information and insight relevant to important applications. For example, in the field of heterogeneous catalysis, the density of states projected onto the d-band at metallic surfaces is a powerful descriptor of a surface's bonding ability, which can be used to estimate its efficiency as a catalyst for a given reaction.
Computational catalysis researchers have found scaling relations between adsorption energies of different adsorbates on metallic surfaces and scaling relations between adsorption energies and reaction barriers that lead to the so-called volcano plots of activity [9]. These volcano plots show that an optimum catalyst should have a 'bonding strength' to atoms and molecules involved in the reaction that enables the reactants to bind to the surface and be activated, but that does not hinder the desorption of the products [10,11]. Concurrently, the d-band model [12] shows a correlation between adsorption energies strengths and some properties of the d-band of a metallic surface, such as its centre [13], width [14], and upper edge [15]. These findings allow using information from the d-projected density of states as indications about how a metallic surface would behave as a catalyst. Such structure-property relationships are pivotal for rational catalyst design, showing the importance of l-p-DOS calculations for this research field.
Our focus in this paper is the ONETEP linear-scaling DFT package, which has recently been used for applications that benefit from l-p-DOS calculations, including studies with metallic nanoparticles applied to heterogeneous catalysis [16][17][18][19]. In the ONETEP linear scaling DFT formulation, the density matrix is represented via non-orthogonal generalised Wannier functions (NGWFs) that are optimised in situ. These atom-centered NGWFs have previously been used to define a l-DOS scheme [20] which has found widespread use. However, after optimisation, the NGWFs are no longer pure in terms of their angular momentum decomposition. Thus, it is necessary to develop and apply a projection scheme to obtain angular momentum resolved properties. Here, we discuss the ONETEP formalism, our choices for the angular momentum resolved basis set, and the implementation of the l-p-DOS method in the ONETEP code. We present test cases, comparing different angular momentum resolved basis sets used to construct the l-p-DOS and comparing the results obtained with ONETEP to similar calculations performed with the CASTEP plane-wave code [21].

The ONETEP code
ONETEP (Order-N electronic total energy package) [22] is a DFT package designed to scale linearly with the number of atoms in the system while retaining some advantages of plane-wave approaches, such as a systematic means to improve the accuracy of the underlying basis. In ONETEP the density matrix is constructed as where φ are the localised non-orthogonal generalised Wannier functions (NGWFs) and K αβ is the density kernel matrix. The Kohn-Sham orbitals are related to the NGWFs by where M α i are the expansion coefficients for the φ α NGWFs. The NGWFs are themselves expanded in terms of a basis of psinc functions [23] which are distributed over the points of a regular grid in real-space whose spacing is controlled by a single parameter. The energy minimisation in ONETEP is performed in two nested loops, one that optimises the density kernel for a fixed set of NGWFs and the other that optimises the NGWFs in situ, allowing the control of the basis set accuracy with a minimal number of NGWFs.
In addition to the locality imposed to the NGWFs, ONETEP also exploits the concept of electronic matter nearsightedness [24] to control the sparsity of the density kernel matrix and pursue linear-scaling computational cost.

Projected density of states in ONETEP
The density of states (DOS) is a function of the energy ( ) that can be defined as: Using the relation from equation (2), we can express this in terms of NGWFs as To construct an angular momentum projected density of states from ONETEP's NGWFs, we first need to define an angular momentum resolved basis, project the eigenfunctions obtained with the NGWFs to this representation and construct a weights matrix that allows summation over particular angular momentum and atoms to decompose the total DOS function. In the following sections, we show the possible choices of angular momentum resolved basis and the methodology used to construct l-p-DOS implemented in ONETEP.

Choice of angular momentum resolved basis
The angular momentum resolved basis functions all have the form of a spherical harmonic, Z l,m (Ω) multiplied by a radial part, where Ω stands for the angular dependence (solid angle). The form of the radial part is where we have some flexibility. Provided that we choose functions which capture enough of the character of the NGWF representation, resulting in a small spilling parameter and a small number of basis functions in the set, the set will be considered adequate. The definition of a 'small' spilling parameter and the number of basis functions considered to be to many or too few is arbitrary and will be explored from now on.
We have explored four options of angular momentum resolved basis, where three of them are based on spherical waves (SW) basis sets. We have as options a non-contracted, full set of SWs; a contracted set with unity contraction weights (C-SWs), a contracted set with contraction weights defined through an inner product with each NGWF, our so-called fitted set (CF-SWs); and finally a set of pseudo-atomic orbitals (PAOs) which is the only set not based directly on SWs, although ultimately the PAOs are themselves expressed via coefficients multiplying the full set of SWs. The pseudo-atomic orbitals are already in use within ONETEP as an starting point from which the NGWFs are optimised, so may represent a low spilling option for the projection of the NGWFs.
We choose the spherical wave (SW) basis as an option so that we can tune the size of the basis systematically to reduce the amount of spilling to the appropriate levels for the application of each calculation. As listed before, this leaves us with some choice, however, in how (or if) we choose to contract the basis functions. Contraction is desirable mainly because our SW basis is not only resolved in the azimuthal quantum number, l, giving the desired angular momentum resolution, but also in magnetic quantum number m and implicitly in the principal quantum number, through the bessel functions, j l . The number of basis functions in this set may be prohibitive for practical calculations, so we can contract over m and j l .
The spherical waves are generated on the same equispaced grid as we use for NGWFs and are defined as: where A represents an atom, H is a Heaviside step function which cuts off any contribution to the spherical waves outside of a radius, a, and r is the vector from the nucleus of the corresponding atom. For each SW, k nl is chosen such that k nl a is the nth zero of the spherical Bessel function j l (x). Contracted spherical waves (C-SWs) are defined in terms of these SWs as: where ω Ak nl lm are the contraction coefficients. For the set denoted C-SWs these are set to 1.
The contraction coefficients can also be set to fit the NGWFs, generating a set of angular momentum resolved functions that we call contracted and fitted spherical waves (CF-SWs). However, in the case that we fit to the NGWFs, we instead take a set of l and m CF-SWs on every NGWF of every atom, rather than one set per atom. This typically increases the size of the spherical wave basis by a factor of 4 for elements up to (and including) the third row of the periodic table and by a factor of 13 for the first transition series, and so on for heavier elements. In this case, the CF-SWs are defined as: where contraction coefficients ω αk nl lm which best fit the NGWF φ α can be calculated by taking the inner product of each NGWF with each SW on the same center: It is possible in this approach that contraction coefficients on an NGWF sum to zero, which will lead to a nonpositive definite overlap matrix, if left untreated. One approach to overcome this problem would be to remove this CF-SW from the set, but instead we opt to leave it in place, recording its index and dealing with it in the resolution of identity, as we explain in the next section.
Another option that we have is to use pseudoatomic orbitals (PAOs) as our angular momentum resolved basis. PAOs are solutions to the atomic Kohn-Sham equation and have the form: where B ln,ν (r) are normalised spherical Bessel functions and c n,ν are a set of coefficients obtained by solving the Kohn-Sham equation for each atom species with spherical localisation constraints using the same cutoff radius as the one used for the NGWFs. Since the PAOs are used for the initialisation of NGWFs [25], to also use them as the angular momentum resolved basis requires no extra calculations associated with the generation of the set. The number of PAOs per atom, N, is the same as the number of NGWFs, typically a number associated with subshell-filling (i.e. 1,4,9,13,...).

Resolution of identity and construction of the projected DOS
Once we have chosen and constructed an angular momentum resolved basis set, we can then construct an identity operator. This can be expressed in terms of the basis functions and their inverse overlap matrix Λ αlm,βl m as: A 'good' angular momentum resolved basis set χ αlm for a given set of NGWFs is implicitly defined by the requirement that in which case it will be safe to insert this resolution of the identity into the density of states from equation (4) without incurring significant 'spillage'. We thus obtain: The major advantage of defining the DOS in terms of this identity operator is that it allows us to exploit the angular momentum dependence of |χ αlm , by selectively summing over particular angular momentum components l and m to create the angular momentum projected DOS.
To calculate the angular momentum projected density of states, we firstly rewrite equation (12) in terms of a weights matrix, W αlm,i : where W αlm,i is defined to be the Hadamard product of two matrices defined as a left and a right part of the identity matrix expression in equation (12), such that and under the condition that The effect of this is to change the order of operations with respect to equation (12), so that instead of producing an element in the identity matrix by multiplying each element of a row in the R T matrix by each element of a column in the T matrix and summing all products in the row and column pair, we instead construct an element in the W matrix by taking the elementwise product of R and T and summing over the elements in the columns of W. This allows for the final sums to be taken over restricted sets of indices, as we require.
For instance, the angular momentum projected density of states can be written as or the local angular momentum projected density of states can be written as where α(S ) means the subset of α associated with the set of atoms, S . Finally, we note that the R αlm i and T αlm,i matrices can be defined in many ways. In our implementation, we choose to define the matrices using a symmetric decomposition achieved by taking the Löwdin decomposition of Λ αlm,βl m , which can be expressed in terms of its eigenpairs, as: where Q is the matrix of eigenvectors of Λ and λ 1 2 is the diagonal matrix of square-rooted eigenvalues of Λ. The T matrix in the symmetric decomposition may then be expressed as where T αlm,i has invariant indices in this case, and W αlm,i from equation (14) can be defined as its elementwise square For all the implemented angular momentum resolved basis, we can define the spilling parameter s, as this helps to quantify how good the angular momentum projected basis is in terms of the accuracy of the identity operator. We will use the spilling parameter later in the paper to assess the quality of our sets. The problem of having basis functions with zero weight leading to eigenvalues of zeros in the C-SWs overlap matrix, discussed in section 3.1 may be dealt with by not inverting Λ αlm,βl m to form Λ αlm,βl m or Λ − 1 2 αlm,βl m , but instead by taking its eigenvalue decomposition which is equivalent to dropping basis functions which lead to non-positive definite overlap matrices. This operation requires a diagonalisation, however, if we need to avoid this, we will need to drop basis functions with zero weights and to use an iterative inversion or Löwdin algorithm. The projected weights allows us to construct various types of DOS. For example, to construct the DOS for the p electrons of a subset of atoms, we can add up all the weights for l = 1 and restrict the summation to the selected atoms: This expression was derived for the case of norm-conserving pseudopotentials where ψ i are the pseudo wavefunctions. In the case of PAW calculations, equation (3) becomes where Ŝ is the PAW overlap operator defined in terms of the PAW transformation operator τ between pseudo and all-electron wavefunctions as Ŝ =τ †τ , and ψ i are the smooth valence wavefunctions.

Calculation details
We performed our simulations with the ONETEP code [22] and with CASTEP [21] as a reference code to validate our implementation. For metallic systems in ONETEP, we used the ensemble DFT method, as implemented by Serrano and Skylaris [26]. For all calculations, we adopted PBE [27] as our exchange-correlation functional. The ONETEP calculations were conducted using the projector augmented wave (PAW) method [28], as implemented in ONETEP by Hine [29], while in CASTEP the calculations were conducted with ultrasoft pseudopotentials [30]. For all the calculations, we have used the GBRV library for ultrasoft and PAW potentials [31]. We have used the OptaDOS [32] code to post-process the results from CASTEP calculations and obtain the l-p-DOS data and set the mid-gap level from ONETEP calculations as zero for all l-p-DOS plots.
In ONETEP, we set the psinc basis set [23] kinetic energy cutoff to 550 eV for geometry optimisations and 850 eV for properties calculations, with the NGWFs radii set to 9.0 a 0 , which has been previously shown to give good convergence for these systems [16][17][18][19]33]. For CASTEP, the kinetic energy cutoff was equal to 450 eV. We performed our simulations under periodic boundary conditions, with a minimum vacuum gap of 10 Å between the borders of the simulation box and any simulated atom. The somewhat lower kinetic energy cutoff required in CASTEP for good convergence compared to ONETEP is largely the result of the need, in ONETEP, for a good description of the discontinuity at the edge of NGWFs spheres. Note also that there is not a direct correspondence between equivalent cutoff values in CASTEP and ONETEP, because the plane-wave basis in CASTEP includes a sphere of G-vectors up to a cutoff defined by 1 2 | G| 2 < E cut , whereas ONETEP uses a real-space grid with a cutoff defined by the same criterion, so implicitly the whole box of G-vectors is used as the basis.
The structures used to compute the projected density of states were optimised in ONETEP [34], where the geometries were relaxed with a convergence threshold of 0.002 Eh/a 0 on the atomic forces. The only exceptions are the Pt slabs and nanoparticles. For these cases, the structures were constructed with bond lengths from the optimised Pt bulk geometry. For all the calculations, the SWs, C-SWs, and CF-SWs were created using a total of 16 functions per magnetic quantum number m per atom.

Demonstration of the method
The spilling parameter presented in equation (22) measures the ability of the basis set to represent the whole manifold of calculated states. Since we aim to compare our basis sets against the equivalent spilling parameter calculated within plane wave methods such as CASTEP, we require a further modification, since plane-wave methods generally restrict the spilling parameter calculation to the occupied bands [7]. A general form of the charge spilling parameter able to deal with fractional occupancies can be written as: where f i are the occupancies and W αlm,i the weights matrix as defined in equation (15).
To compare the quality of the angular momentum resolved bases, we present in table 1 the charge spilling parameter for different systems. We calculated the l-p-DOS for a CO molecule, a set of hydrocarbons, a silane molecule, a LiPF 6 compound in the C 4v symmetry, a set of cuboctahedral metallic nanoparticles (Pt and Pd) and a platinum slab with 320 atoms and a (1 1 1) facet exposed. We performed the calculations using all the available l-p-DOS bases in ONETEP, namely: (i) non-contracted spherical waves (SWs); (ii) contracted spherical waves (C-SWs); (iii) contracted spherical waves with coefficients fitted to the NGWFs (CF-SWs) and; (iv) pseudoatomic orbitals (PAOs). Table 1 shows a general trend in the charge spilling parameters, namely that SWs < PAOs CASTEP PAOs < CF-SWs < C-SWs. The small values obtained with SWs were expected due to the large size and flexibility of the basis. The results obtained with PAOs in ONETEP and CASTEP indicate that the implementation in both codes is comparable and that we should also expect similar results for pDOS calculations. We can also note that the spilling parameters from CASTEP are slightly higher than the ones obtained with PAOs in ONETEP, which can be seen as a consequence of the localisation of the NGWFs in spheres with the same radius as the PAOs, whereas in CASTEP the wavefunctions are delocalised over the simulation box. Table 1 also shows that determining the contraction coefficients by fitting the C-SWs to the NGWFs contributes to reducing the spilling parameter as compared with contraction of spherical waves using unit weights. For the calculations with metallic nanoparticles, all the pDOS options, with the exception of C-SWs, show similar values of spilling parameters, indicating that the projected density of states obtained with these approaches should also be comparable to each other.
We compare in figure 1 projected density of states obtained with CASTEP and the ones obtained with PAOs in ONETEP. Figures 1(a) and (b) show, respectively, the density of states projected onto s and p orbitals for a C 2 H 4 molecule. Figure 1(c) tests, simultaneously, the local and angular momentum projected density of states, by showing the DOS projected onto the d-channel only for Pt atoms in a (1 1 1) facet of a cuboctahedral Pt 13 nanoparticle.
In accordance with the results obtained for the charge spilling parameters, here we see a very high degree of agreement between ONETEP and CASTEP results. This is also observed for several other tested systems but omit- Table 1. Charge spilling parameter, presented as percentages, for calculations performed with CASTEP, which uses its own set of PAOs as the basis set for of its charge spilling calculations, and different angular momentum resolved bases in ONETEP. From now on, we will use only the results obtained with the PAOs in ONETEP to compare with other basis sets. Figure 2 compares all the implemented angular momentum resolved basis sets, testing the results for CO, C 2 H 4 and C 2 H 6 . Figures 2(a), (c) and (e) show, respectively, the density of states of CO, C 2 H 4 and C 2 H 6 projected onto s channels, while (b), (d) and (f) shows the projections for p channels of the same molecules. In general, the projections with all the angular momentum resolved bases are similar. The generation of the non-contracted spherical waves set is unbiased and, as demonstrated with the charge spilling parameters, provides an almost complete projection of the calculated bands for a wide range of systems. In this sense, the similarity of the pDOS obtained with PAOs and non-contracted SWs reinforces the idea that PAOs are viable options to obtain angular momentum projected density of states.
The main differences arise when calculating systems with atoms with valence shells composed by different values of azimuthal quantum numbers. The spherical wave basis is constructed with the same maximum l for every atom in the whole system. For example, for a hydrocarbon, the SWs are constructed with s and p functions for carbon and hydrogen atoms. Meanwhile, the PAOs are constructed by default in ONETEP as a minimum Out of all the contracted spherical wave schemes, the CF-SWs is the one that presents largest variations from PAOs and non-contracted SWs. We attribute this to the difficulty of fitting a spherical wave basis composed of s and p functions to a single NGWF. Thus, we have tested the effect of including five NGWFs on hydrogen atoms, so as to initialise them with s (2 NGWFs) and p (3 NGWFs) character. As a result, the calculated charge spilling parameters decrease from 7.44% to 3.69% for C 2 H 4 and from 6.96% to 2.33% for C 2 H 6 , and the overall behaviour observed with non-contracted spherical waves is recovered.
As briefly discussed in the introduction, angular momentum projected density of states is an invaluable tool for theoretical investigations in the field of heterogeneous catalysis. In this area, describing the interaction between the catalyst surface and adsorbates is essential to explain and predict catalytic activities. As proposed by Nørskov et al [12,13], properties of the localised d-band of transition metals can be used as electronic descriptors of the interaction strength between and adsorbate and the catalyst surface. Specifically, as there is a connection between the d-band centre and the surface composition, crystallographic orientation, site coordination and nanoparticle sizes, studies with such descriptors can be used to search for directions to catalysts optimisations.
Here, to test the suitability of our l-p-DOS methods for application to this area, we have decided to compute cuboctahedral platinum nanoparticles of increasing size, ranging from Pt 13 to Pt 561 and a Pt slab with a (1 1 1) facet exposed. We have also computed a Pd 13 cluster as an additional test. Increasing the nanoparticle size affects the Pt-Pt bond lengths at the surface, and the ratio between uncoordinated and coordinate sites. This generates shifts in the d-band centres that translate into a weakening in the interaction with adsorbates. Here, we are not considering the changes in the Pt-Pt bond lengths with increasing nanoparticle sizes. We compute Pt nanoparticles created with Pt bulk bond lengths to isolate the electronic effects to the geometrical changes due to size effects. Figures 3(a) and (b) present the density of states projected onto the d channels of Pt atoms in a (1 1 1) facet of a Pt 147 nanoparticle. The first plot spans all the calculated energy levels and the second focuses on the states near the Fermi level. These plots show that d-bands calculated with all the implemented basis sets are similar. The overall shape of the d-band near the Fermi level remains almost unchanged for all basis sets, with only small variations in terms of peak intensities. The main differences appear for low energy levels, where additional d-projected states are observed for spherical wave based approaches. In practice, as the energy levels near the Fermi level are the important region for catalysis applications, and these results are similar for all basis set options, these differences should not affect the overall applicability of any of the schemes. Figure 4(a) shows the d-band centres calculated with different approaches for the (1 1 1) facet of the nanoparticles against the nanoparticles sizes. The d-band centres are calculated using the occupancy-weighted l-p-DOS and focusing on the states near the Fermi-level, i.e. excluding the variations for low energy levels commented on figure 3. As we increase the nanoparticle size, the d-band centre moves away from the Fermi-level, which according to the d-band theory should describe a weaker interaction with adsorbates. Figure 4(b) shows how the calculated d-band centres with spherical waves bases (SWs, C-SWs, CF-SWs) compare to the same quantities calculated using PAOs. We show the d-band centres for the whole (1 1 1) facet of Pt nanoparticles of increasing size, the exposed facet of the Pt (1 1 1) slab and the (1 1 1) facet of Pd 13 . As observed, the small changes in the projected density of states plots are not crucial for the calculated d-band centres and similar trends and conclusions for the catalytic activity of the nanoparticles could be obtained with any method.

Conclusions
In this work, we have described the underlying theoretical framework for the implementation of angular momentum projected density of states within the ONETEP formalism. We presented four options of angular momentum resolved basis sets that were implemented in ONETEP to obtain l-p-DOS. The first three are based on spherical waves, with a non-contracted set of SWs; a contracted set with unity contraction weights (C-SWs), and a contracted set where the contraction weights are determined via inner products with the NGWFs (CF-SWs). Pseudo-atomic orbitals (PAOs) are used for the fourth basis set, which is also the initialisation of NGWFs in ONETEP.
We performed tests on several systems to assess how the projection varies with changes in the basis sets and how the ONETEP's results compare to l-p-DOS functionality in CASTEP. We have calculated the density of states for a CO and a silane molecule, a set of hydrocarbons, a LiPF 6 compound, metallic nanoparticles of different sizes, and one platinum slab with 320 Pt atoms. In terms of charge spilling parameter, we observed a similar trend for all tested systems, with SWs < PAOs CASTEP < CF-SWs < C-SWs. For the non-contracted SWs basis, small charge spilling parameters s q < 1% were observed for all systems. For PAOs, the charge spilling parameters obtained with ONETEP and CASTEP were almost identical, and for contracted SWs basis, the option with contraction weights obtained by fitting the NGWFs presented lower s q values than the one with unit contraction weights.
Despite some differences in the charge-spilling parameters, the l-p-DOS results are comparable between the implemented approaches and between the two tested codes. The more evident differences arise with systems composed of atoms with different values of azimuthal quantum numbers in the valence shell, where the l-p-DOS obtained with CF-SWs diverge from the other options. This difference can be reduced by adding NGWFs in the ONETEP basis set to allow a better fitting between NGWFs and SWs. Moreover, for metallic nanoparticles of increasing size, the d-bands near the Fermi level, which is commonly used as a descriptor in the heterogeneous catalysis field, remain almost unaltered between different approaches, showing the robustness of the implemented method for such studies. The availability of these methods in a linear-scaling framework such as ONETEP opens the way for analysis of the electronic structure of complex nanostructured materials.