The Wannier-projected density of states, populations and charges

As well as enabling the calculation of bonding descriptors from Wannier functions, pengwann also allows for the computation of a Wannier-projected density of states which can be integrated to derive atomic populations and charges. This example demonstrates the calculation of these quantities in rutile TiO\(_{2}\).

The Wannier-projected density of states

from pengwann.geometry import Geometry

geometry = Geometry.from_xyz(seedname='wannier90', path='inputs')

print(geometry)

Geometry
========
Cell
----
[[4.587891 0.       0.      ]
 [0.       4.587891 0.      ]
 [0.       0.       2.938132]]

Assignments
-----------
X0 => (28,)
X1 => (28,)
X2 => (28,)
X3 => (28,)
X4 => (28,)
X5 => (28,)
X6 => (29,)
X7 => (29,)
X8 => (29,)
X9 => (29,)
X10 => (29,)
X11 => (29,)
X12 => (30,)
X13 => (30,)
X14 => (30,)
X15 => (30,)
X16 => (31,)
X17 => (31,)
X18 => (31,)
X19 => (31,)
X20 => (32,)
X21 => (32,)
X22 => (32,)
X23 => (32,)
X24 => (33,)
X25 => (33,)
X26 => (33,)
X27 => (33,)
Ti28 <= (0, 1, 2, 3, 4, 5)
Ti29 <= (6, 7, 8, 9, 10, 11)
O30 <= (12, 13, 14, 15)
O31 <= (16, 17, 18, 19)
O32 <= (20, 21, 22, 23)
O33 <= (24, 25, 26, 27)

In previous examples, we have used the identify_interatomic_interactions() function to find bonds and describe them in terms of the interacting atoms’ assigned Wannier functions. To instead calculate local atomic properties such as the Wannier pDOS, we make use of the identify_onsite_interactions() function:

from pengwann.geometry import identify_onsite_interactions

# The atomic species for which to obtain interactions.
symbols = ('Ti', 'O')

interactions = identify_onsite_interactions(geometry, symbols)

print(interactions)

Atomic interactions
===================
Ti28 <=> Ti28
Ti29 <=> Ti29
O30 <=> O30
O31 <=> O31
O32 <=> O32
O33 <=> O33

The meaning of “on-site” interactions in this context is made slightly clearer if we examine a specific AtomicInteraction:

print(interactions[28, 28])

Atomic interaction Ti28 <=> Ti28
================================
DOS matrix => Not calculated
WOHP => Not calculated
WOBI => Not calculated
IWOHP => Not calculated
IWOBI => Not calculated
Population => Not calculated
Charge => Not calculated


Associated Wannier interactions
-------------------------------
0[0, 0, 0] <=> 0[0, 0, 0]
1[0, 0, 0] <=> 1[0, 0, 0]
2[0, 0, 0] <=> 2[0, 0, 0]
3[0, 0, 0] <=> 3[0, 0, 0]
4[0, 0, 0] <=> 4[0, 0, 0]
5[0, 0, 0] <=> 5[0, 0, 0]

As seen above, an on-site interaction in pengwann refers to all of the diagonal terms, in which the coefficients used to calculate the DOS matrix are related by complex conjugation, yielding the Wannier-projected density of states:

\[\mathrm{pDOS}_{\alpha}(E) = \sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\alpha}_{nk}\right]\cdot\delta(E - \epsilon_{nk}).\]

Note

For more details regarding the mathematical formalism employed in pengwann, please see the methodology page and the API reference.

from pengwann.descriptors import DescriptorCalculator
from pengwann.io import read

kpoints, eigenvalues, u_matrices, _ = read(seedname='wannier90',
                                           path='inputs')

mu = 9.5
nspin = 2

num_wann = 28
energy_range = (-20, 20)
resolution = 0.1
sigma = 0.1

dcalc = DescriptorCalculator.from_eigenvalues(eigenvalues,
                                              num_wann,
                                              nspin,
                                              energy_range,
                                              resolution,
                                              sigma,
                                              kpoints,
                                              u_matrices)

Having initialised a DescriptorCalculator in much the same way as we would to calculate bonding descriptors, we can now obtain the pDOS by calling the assign_descriptors() method on our on-site interactions:

interactions = dcalc.assign_descriptors(interactions,
                                        calc_wohp=False,
                                        calc_wobi=False)

Note that, as we are only interested in the DOS matrix in this case, we have set the calc_wohp and calc_wobi keyword arguments to False to avoid unnecessarily fetching or computing elements of the Wannier Hamiltonian and Wannier density matrix.

We can now plot the pDOS for a given atom simply by accessing the dos_matrix attribute. For example, for Ti28:

import matplotlib.pyplot as plt

# Colours for plots
ti_colour = '#79a9d1'
o_colour = '#ff4365'

shifted_energies = dcalc.energies - mu

fig, ax = plt.subplots()

pdos = interactions[28, 28].dos_matrix
ax.plot(shifted_energies, pdos, color=ti_colour)
ax.axvline(x=0, color='black', ls='--', lw=1)

ax.set_xlim(-20, 8)

ax.set_xlabel(r'$E - E_{\mathrm{F}}$ / eV')
ax.set_ylabel('pDOS')

plt.show()

Note

In much the same way that bonding descriptors such as the Wannier orbital Hamilton population (WOHP) can be resolved with respect to individual Wannier functions (see this example), the same is true of the pDOS.

To plot the pDOS in terms of atomic species rather than individual atoms, we simply sum over atoms of the same species. Accessing all interactions involving a certain subset of atomic species/elemental symbols is easily accomplished using the filter_by_species() method:

ti_pdos = sum(
    interaction.dos_matrix for interaction in 
    interactions.filter_by_species(['Ti'])
)
o_pdos = sum(
    interaction.dos_matrix for interaction in 
    interactions.filter_by_species(['O'])
)

fig, ax = plt.subplots()

ax.plot(shifted_energies, ti_pdos, color=ti_colour, label='Ti')
ax.plot(shifted_energies, o_pdos, color=o_colour, label='O')
ax.axvline(x=0, color='black', ls='--', lw=1)

ax.set_xlim(-20, 8)

ax.set_xlabel(r'$E - E_{\mathrm{F}}$ / eV')
ax.set_ylabel('pDOS')

ax.legend(frameon=False)

plt.show()

Populations and charges

The pDOS can be integrated to derive a Löwdin-style population for a particular Wannier function or atom.

interactions = interactions.with_integrals(dcalc.energies, mu)

for interaction in interactions:
    print(f'Population for {interaction.tag} = {interaction.population:.2f}')

Population for Ti28 <=> Ti28 = 1.88
Population for Ti29 <=> Ti29 = 1.88
Population for O30 <=> O30 = 7.06
Population for O31 <=> O31 = 7.06
Population for O32 <=> O32 = 7.06
Population for O33 <=> O33 = 7.06

If we also provide valence_counts to the with_integrals() method, then we will also have access to charges for each atom too:

valence_counts = {'Ti': 4, 'O': 6}

interactions = interactions.with_integrals(dcalc.energies,
                                           mu,
                                           valence_counts=valence_counts)

for interaction in interactions:
    print(f'Charge for {interaction.tag} = {interaction.charge:.2f}')

Charge for Ti28 <=> Ti28 = 2.12
Charge for Ti29 <=> Ti29 = 2.12
Charge for O30 <=> O30 = -1.06
Charge for O31 <=> O31 = -1.06
Charge for O32 <=> O32 = -1.06
Charge for O33 <=> O33 = -1.06

The relationship between populations and charges in this instance is trivial, the charge is simply the difference between the number of valence electrons associated with an atom (according to the pseudopotential used in the ab initio calculation) and its population.