# Calculating the Local/Partial Density of States and Angular Momentum Projected Density of States

- Author:
Nicholas D.M. Hine, University of Warwick (originally Imperial College London)

- Author:
Jolyon Aarons, University of Warwick

- Author:
Weibo Ng, University of Southampton

- Date:
August 2024 (Updated by Weibo Ng to include recommended parameters).

- Date:
June 2019 (Updated by Jolyon Aarons to add angular momentum PDOS information).

- Date:
Originally written by Nicholas D.M. Hine April 2012.

## What is being calculated?

The local density of states (LDOS) provides a description encompassing both the spatial and energetic distribution of the single-particle eigenstates simultaneously. The angular momentum projected density of states (PDOS) decomposes the density of states energetic distribution into angular momentum components. Both decompositions may be combined into an LPDOS. The LDOS and PDOS can thus be two valuable sources of information for understanding and interpreting electronic structure calculations.

In the local-orbital framework of ONETEP [Skylaris2005], both decompositions are achieved by first performing a diagonalisation of the Hamiltonian matrix in the NGWF basis. This is a post-processing step performed at the end of the calculation, once NGWF and density kernel convergence have been achieved. While this comes with a \(O(N^{3})\) computational cost, the prefactor is low because the NGWF basis is generally quite small. Therefore, such a diagonalisation remains fast up to quite large system sizes, particularly if a parallel eigensolver such as ScaLAPACK is used.

The generalised eigenproblem that needs to be solved to provide the eigenvalues and eigenvectors is:

The matrix \(M_{\phantom{\beta}n}^{\beta}\) describes the eigenvectors, which take the form \(|\psi_{n}\rangle=\sum_{\beta}|\phi_{\beta}\rangle M_{\phantom{\beta}n}^{\beta}\). In ONETEP, Eq. (52) can be solved using either the LAPACK routine DSYGVX or (preferably) the ScaLAPACK routine PDSYGVX, depending on whether -DSCALAPACK has been specified at compile time, and ScaLAPACK libraries have been provided.

The result is the eigenvalues \(\{\epsilon_{n}\}\) and eigenvectors \(M_{\phantom{\beta}n}^{\beta}\). From these, the total density of states can be obtained as

In practice the delta function is replaced with a Gaussian with broadening \(\sigma\), typically of the order of \(0.1\,\mathrm{eV}\):

### Local Density of States

The local density of states in a given region \(I\) is calculated by projecting each eigenstate onto the local orbitals of region \(I\), as

Here, the non-orthogonality of the NGWFs when used as projectors means that the ket must be contravariant. We are therefore implicitly using the contravariant dual of the NGWF (as in DFT+U [O-Regan2011], [O-Regan2012]). Fortunately, the functions \(|\phi^{\alpha}\rangle\) need not be explicitly constructed in real space: the relationship \(\langle\phi_{\alpha}|\phi^{\beta}\rangle=\delta_{\alpha\beta}\) implies that we can re-write Eq. (54) as:

What is therefore obtained is a series of functions \(D_{I}(\epsilon)\) for each of the chosen regions \(I\), which may be the NGWFs of a single atom, or those of a group of atom types.

### Angular Momentum Projected Density of States

To calculate the PDOS, we need to insert an additional resolution of the identity into equation (55) using a basis of angular momentum resolved functions on which to project our NGWFS, \(| \chi_{l,m}\rangle\):

where we need to include the overlap matrix of angular momentum resolved functions, \(\Lambda\), since this basis is also non-orthogonal.

We have considerable scope in which basis we choose for the angular momentum resolved functions. Effectively, this is a set of spherical harmonics multiplied by some radial term. In ONETEP, we currently have two options implemented for the radial term: either spherical waves, or pseudo-atomic functions, as used to initialise the NGWFs, before optimisation in the NGWF SCF loop. More details about the theory behind these options as well as tests and comparisons can be found in our paper [Aarons2019].

## Performing an LDOS Calculation

An LDOS calculation is performed as part of the optional post-processing
activated using `do_properties: T`

or using `task: PROPERTIES`

.
To activate LDOS we then need to specify the Gaussian broadening, such
as `dos_smear : 0.1 eV`

. The default value of
`dos_smear : -0.1 eV`

disables LDOS.

Then, we need to specify the groups of atom types. This is done via a block, with each line listing a group of atoms. For example, in a benzene ring, we might use the following to find the contributions of the carbon and hydrogen atoms respectively:

```
%block species_ldos_groups
C
H
%endblock species_ldos_groups
```

A more complex example would be for a GaAs nanorod with hydrogen termination on the faces. If we wished to see the LDOS varying over 5 layers, labelled 1-5, we could use:

```
%block species_ldos_groups
Ga1 As1 H1
Ga2 As2 H2
Ga3 As3 H3
Ga4 As4 H4
Ga5 As5 H5
%endblock species_ldos_groups
```

Examples of the use of LDOS analysis, including example plots, can be found in several recent papers employing ONETEP [Avraam2011], [Avraam2012], [Hine2012].

## Performing a PDOS Calculation

The default settings in ONETEP for PDOS calculations are to use the
pseudo-atomic states as the angular momentum resolved projection basis
with a Löwdin orthogonalisation. For most applications, the spilling
parameter associated with this basis will be sufficiently small.
PDOS calculations are enabled as part of the optional post-processing by
writing `do_properties : T`

into the ONETEP input file, along with a
Gaussian smearing width, such as `dos_smear : 0.1 eV`

and a maximum
angular momentum in the angular momentum resolved projection basis, such
as `pdos_max_l : 2`

to include up to d-states (when using the
default pseudo-atomic basis, this will also be limited by the maximum
angular momentum state in each species, as calculated by the
pseudo-atomic solver).

### Local PDOS (LPDOS) Calculations

If you intend to calculate an angular momentum projected DOS on a subset
of atoms, this can be achieved by specifying a block in the input file,
in the same way as for LDOS. The block can be set up by using the
`species_pdos_groups`

keyword. For example, in a benzene ring
calculation, if you want to find the contributions to the PDOS coming
from solely carbon atoms and solely hydrogen atoms, you could write:

```
%block species_pdos_groups
C
H
%endblock species_pdos_groups
```

If you also want the combined contribution from carbon and hydrogen to each angular momentum channel, you should add a line for this:

```
%block species_pdos_groups
C
H
C H
%endblock species_pdos_groups
```

This will calculate PDOS histogram data up to `pdos_max_l`

for each
line. As many or as few combinations of species as you require can be
calculated by adding extra lines.

If you instead want a specific subset of atoms of a particular species, this can be achieved easily by labelling this subset differently to the others in its species. For example, if you have the following species specification:

```
%block species
Pt Pt 78 9 9.0
%endblock species
%block species_atomic_set
Pt SOLVE conf=5d9 6s1 6p0
%endblock species_atomic_set
%block species_pot
Pt platinum.paw
%endblock species_pot
%block species_pdos_groups
Pt
%endblock species_pdos_groups
%block positions_abs
Pt 8.8292 12.2847 8.7330
Pt 9.2819 11.1839 11.1325
%endblock positions_abs
```

`positions_abs`

block, as shown
here for example:```
%block species
Pt Pt 78 9 9.0
Pt1 Pt 78 9 9.0
%endblock species
%block species_atomic_set
Pt SOLVE conf=5d9 6s1 6p0
Pt1 SOLVE conf=5d9 6s1 6p0
%endblock species_atomic_set
%block species_pot
Pt platinum.paw
Pt1 platinum.paw
%endblock species_pot
%block species_pdos_groups
Pt
Pt1
Pt Pt1
%endblock species_pdos_groups
%block positions_abs
Pt 8.8292 12.2847 8.7330
Pt1 9.2819 11.1839 11.1325
%endblock positions_abs
```

### Expert PDOS Options

Further options are available in the ONETEP PDOS functionality to control the quality of the projection. These will be unneeded in most cases, but if, for instance, you are observing larger spilling parameters than your requirements permit, you may wish to enable some of these options.

The most reliable way we have found to reduce the spilling parameter is
to use a spherical-wave basis rather than the pseudo-atomic basis as the
angular momentum resolved projection basis. To do this in ONETEP, add
`pdos_pseudoatomic : F`

to your input file. By default, this will
create a set of contracted spherical waves by fitting spherical waves to
your converged NGWFs, via the contraction coefficients.

The spherical wave basis is contracted by default to reduce the memory
requirements of the code. You may, however, not see an improvement in
the spilling parameter by using this set. To be certain of reducing the
spilling parameter, you should also opt to use the full, non-contracted
spherical wave basis, by setting `pdos_reduce_sws F`

in your input
file, along with an adequately large `pdos_max_l`

. For `pdos_max_l`

you can start by running with 2 and increase to 3 if required. If you
choose to take this approach, beware of the memory requirements, which
can be *up to* 10 times greater.

If you choose to use a contracted set, then you almost certainly want to
use the default fitting coefficients (fitted to NGWFs). These can be
changed to unity by setting `pdos_construct_basis : F`

, however this
is not likely to improve your results, and is likely to be removed in a
future version of ONETEP due to it being mainly of use for debugging
purposes. To reduce the spilling parameter with the contracted set, we
recommend increasing the `pdos_max_l`

parameter.

In specialised cases, you may also wish to *not* sum over the magnetic
quantum number. This can be achieved by setting `pdos_sum_mag : F`

.
This will give histogram data for every magnetic quantum number of every
angular momentum channel of each atom group.

### Recommended Settings

The following table provides a set of recommended settings for calculating the PDOS in ONETEP:

Accuracy |
Cutoff energy (eV) |
NGWF radii (\(a_{0}\)) |
Basis |
---|---|---|---|

Low |
700 |
8.0 |
Pseudo-atomic orbitals |

Medium |
800 |
10.0 |
Pseudo-atomic orbitals |

High |
1200 |
13.0 |
Fully uncontracted spherical waves |

These settings were calibrated and tested with platinum nanoparticles and can be considered a starting point for your own calculations.

### Interpreting Outputs

ONETEPs PDOS outputs are written to several files as well as to stdout,
which is itself usually redirected to the main log/output file. The
PDOS output to stdout will look something like the following (for an
input file with 3 `pdos_groups`

and `pdos_max_l=2`

):

```
================ Projected Density of States (pDOS) calculation ================
Constructing AM resolved functions ... done
Writing NGWF plot files in formats: ... done
Performing overlap integrals ... done
Computing pDOS weights ... done
All bands spilling parameter = 2.14 %
Occupancy-weighted spilling parameter = 0.30 %
=> Outputting data for OptaDOS <=
Writing pDOS weights to file "Pt3O.val_pdos_bin" ... done
Writing band gradients to file "Pt3O.val_dome_bin" ... done
Writing Castep output cell file to "Pt3O-out.cell" ... done
=> Computing Gaussian smeared pDOS <=
Writing "Pt3O_PDOS.txt" ... done
=> Computing Occupancy-weighted Gaussian smeared pDOS <=
Writing "Pt3O_occ_PDOS.txt" ... done
=> Band centres:
S band centre of group 1 from -24.409 eV upwards: -10.785237 eV
P band centre of group 1 from -24.409 eV upwards: -6.380437 eV
D band centre of group 1 from -25.000 eV upwards: -1.992329 eV
S band centre of group 2 from -24.409 eV upwards: -3.492322 eV
P band centre of group 2 from -24.409 eV upwards: -5.494682 eV
D band centre of group 2 from -25.000 eV upwards: -1.992329 eV
S band centre of group 3 from -24.409 eV upwards: -20.033099 eV
P band centre of group 3 from -24.409 eV upwards: -6.607392 eV
Band centres done. <=
=> Integrated number of electrons in each AM band:
S num electrons of group 1 from -24.409 eV upwards: 3.768881
P num electrons of group 1 from -24.409 eV upwards: 5.624424
D num electrons of group 1 from -25.000 eV upwards: 26.497489
S num electrons of group 2 from -24.409 eV upwards: 2.107161
P num electrons of group 2 from -24.409 eV upwards: 1.147191
D num electrons of group 2 from -25.000 eV upwards: 26.497489
S num electrons of group 3 from -24.409 eV upwards: 1.661719
P num electrons of group 3 from -24.409 eV upwards: 4.477233
Integrated number of electrons done. <=
================================================================================
```

First, we can see the spilling parameters – effectively how well the angular momentum resolved basis is able to represent the NGWFs. A lower value is better; if you are running production calculations, you should want a value lower than a few percent. If not, consider making some of the changes suggested in the expert options section above. The all-bands value includes un-occupied bands as well as valence states, which will be the same as the occupancy-weighted version unless you are using EDFT with a finite electronic temperature.

Following this are the files for use with **OptaDOS**. OptaDOS is a
freely available piece of software for plotting various DOS projections.
If you intend to use it, then please follow the CASTEP section of the
OptaDOS manual with these files, as they are compatible.

The histogram files are then written, “*_PDOS.txt”. These also come in all-bands and occupancy weighted flavours (the occupancy weighted variant is more reliable for a usual ground state calculation with ONETEP as the conduction states are not well described. Only if you are doing LPDOS on the output of a ONETEP conduction calculation the occupancy un-weighted LPDOS outputs will be meaningful). The order of columns is firstly the energy column, followed by the angular momentum component columns (i.e. s,p,d…) for each pdos group. This can be plotted trivially with xmgrace, or any other plotting tool.

Finally ONETEP reports the energy and occupancy weighted averages of the
PDOS, so called-band centres, useful in catalysis (e.g. the value
“d-band centre” is a very useful descriptor about the ability of a metal
surface to bind atomic oxygen and other types of adsorbates) and the
integrated number of electrons in each component. The d-band centre and
number of electrons in each band are calculated above a threshold.
By default, this is set to -15 eV, which works well for platinum nanoparticles,
but this can be adjusted using the `pdos_d_band_threshold`

keyword.
The threshold was introduced as d-orbital contributions occured at regions
outside of the d-band, which can be seen in Figure 3a. of [Aarons2019].

[Skylaris2005] C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, J. Chem. Phys. **122**, 084119 (2005).

[O-Regan2011] D. D. O’Regan, M. C. Payne and A. A. Mostofi, Phys. Rev. B **83**, 245124 (2011).

[O-Regan2012] D. D. O’Regan, N. D. M. Hine, M. C. Payne and A. A. Mostofi, Phys. Rev. B **85**, 085107 (2012).

[Avraam2011] P. W. Avraam, N. D. M. Hine, P. Tangney, and P. D. Haynes, Phys. Rev. B **83**, 241402(R) (2011).

[Avraam2012] P. W. Avraam, N. D. M. Hine, P. Tangney, and P. D. Haynes, Phys. Rev. B **85**, 115404 (2012).

[Hine2012] N. D. M. Hine, P. W. Avraam, P. Tangney, and P. D. Haynes, J. Phys. Conf. Ser. (2012).

[Aarons2019] J. A. Aarons, L. G. Verga, N. D. M. Hine, and C.-K. Skylaris, Submitted (2019).