Advanced analysis of the electronic structure

In the basic tutorials we have seen ways to “see” the electronic structure via DOS plots, band structures, or wavefunction pictures. Here we take a deeper view, using as basic ingredients the coefficients of the wavefunctions and the matrix elements of the hamiltonian and overlap matrix.

In Siesta, the wave-functions are expanded in localized orbitals:

\[|\Psi_i> = \sum_{\mu} { c^i_{\mu} |\mu>}\]

The coefficients \(c^i_{\mu}\) are produced by the diagonalization routines in Siesta and stored in the *.WFSX files.

The electrons will fill up the states up to the Fermi level. If we denote the occupation of a state i by \(f_i\), we can write for the total band-structure energy:

\[E_{BS} = \sum_i { f_i <\Psi_i|H|\Psi_i>} = \sum_i { f_i \sum_{\mu\nu} { c^i_{\mu}c^i_{\nu} H_{\mu\nu} } }\]

where \(H_{\mu\nu} = <\mu|H|\nu>\) are the matrix elements of the Hamiltonian in the basis set. Alternatively:

\[E_{BS} = \sum_i { f_i <\Psi_i|H|\Psi_i>} = \sum_i { f_i \sum_{\mu\nu} { \varepsilon_i c^i_{\mu}c^i_{\nu} S_{\mu\nu}}}\]

where \(S_{\mu\nu} = <\mu|\nu>\) are the elements of the overlap matrix, which is needed in the case of non-orthogonal orbitals.

In the same way, the total number of electrons is given formally by:

(1)\[N = \sum_i { f_i <\Psi_i|\Psi_i>} = \sum_i { f_i \sum_{\mu\nu} { c^i_{\mu}c^i_{\nu} S_{\mu\nu}}}\]

Mulliken charges

The last equation can be re-arranged as:

\[N = \sum_I { \sum_{\mu\in I} { \sum_{\nu} { \sum_i { f_i c^i_{\mu}c^i_{\nu} S_{\mu\nu}}}}}\]

where we have partitioned the \(\mu\) index in different subsets for different atoms I. It is now a small logical step to think of

\[q_I = \sum_{\mu\in I} { \sum_{\nu} { \sum_i { f_i c^i_{\mu}c^i_{\nu} S_{\mu\nu}}}}\]

as the “contribution of atom I” to the total charge of the system. These are the “Mulliken populations”, and are printed by Siesta if the option Write-Mulliken-Pop 1 is used.

For example, for the water molecule with a DZP basis, we get in the output file:

mulliken: Atomic and Orbital Populations:

Species: O
Atom  Qatom  Qorb
               2s      2s      2py     2pz     2px     2py     2pz     2px
   1  5.716   0.732   0.619   1.311   1.714   0.588   0.173   0.064   0.494

Species: H
Atom  Qatom  Qorb
               1s      1s      2Ppy    2Ppz    2Ppx
   2  1.142   0.392   0.495   0.109   0.110   0.036
   3  1.142   0.392   0.495   0.109   0.110   0.036

mulliken: Qtot =        8.000

which seems to show that the oxygen atom has lost a bit of charge (?), and the hydrogen atoms have gained electrons (??). The problem, quoting the Wikipedia article is:

[…] Mulliken charges are explicitly sensitive to the basis set choice. In principle, a complete basis set for a molecule can be spanned by placing a large set of functions on a single atom. In the Mulliken scheme, all the electrons would then be assigned to this atom. The method thus has no complete basis set limit, as the exact value depends on the way the limit is approached. This also means that the charges are ill defined, as there is no exact answer. As a result, the basis set convergence of the charges does not exist, and different basis set families may yield drastically different results.

In fact, if we use a SZ basis, on paper a much worse one, we get:

Species: O
Atom  Qatom  Qorb
               2s      2py     2pz     2px
   1  6.616   1.803   1.624   2.000   1.189

Species: H
Atom  Qatom  Qorb
               1s
   2  0.692   0.692
   3  0.692   0.692

which is closer to what one expects.

Other populations: Hirshfeld, Voronoi, (Bader)

The simple recipe for the Mulliken populations was too good to be true. There are other ways to tackle the (essentially ill-defined) problem of assigning charge to “atoms” (in the DFT nindset, what we really have is an “electron soup” that has found its optimal arrangement in a field of positively charged ions. One can see “atomic-like” features but well-defined partitions do not exist).

In Siesta we can use the options (…) to request the computation of Hirshfeld and Voronoi charges. The first are based on a partition that uses the footprint of a superposition of free-atom charges, and the second are associated to the Voronoi polyhedron around an atom.

For the first case above, these charges come out to be:

Hirshfeld Atomic Populations:
Atom #   dQatom  Atom pop  Species
   1   -0.23391   6.23391  O
   2    0.11696   0.88304  H
   3    0.11696   0.88304  H

Voronoi Atomic Populations:
Atom #   dQatom  Atom pop  Species
   1   -0.17545   6.17545  O
   2    0.08773   0.91227  H
   3    0.08773   0.91227  H

and qualitatively the same for the SZ case. You can check that these charges are indeed quite insensitive to the details of the basis set.

Finally, Siesta can generate output to be used by analysis programs based in Bader’s “atoms in molecules” scheme. Please consult the manual for more details.

DOS and projected DOS

Now we can take a seemingly crazy step and write the occupation numbers in the above equations as:

(2)\[f_i = \int_0^{\varepsilon_F} { \delta(\varepsilon-\varepsilon_i) d\varepsilon }\]

since the integral will be one for states below the fermi level \(\varepsilon_F\) and zero for those above it. Then, the decomposition of N reads:

\[N = \int_0^{\varepsilon_F} { d\varepsilon \sum_i { \sum_{\mu} { \sum_{\nu} { c_{i\mu}c_{i\nu} S_{\mu\nu} \delta(\varepsilon - \varepsilon_i) } } } }\]

This seemingly pointless way to write N is actually quite useful. Begin by taking out the integral sign:

\[g(\varepsilon) = \sum_i { \sum_{\mu} { \sum_{\nu} { c_{i\mu}c_{i\nu} S_{\mu\nu} \delta(\varepsilon - \varepsilon_i) } } }\]

This is the density of states, which can also be written more simply as

\[g(\varepsilon) = \sum_i { \delta(\varepsilon - \varepsilon_i) }\]

since the sum over \(\mu\nu\) is one for each i (exercise to the reader).

In the next stage, we take out the sum over \(\mu\), and get:

\[g_{\mu} (\varepsilon) = \sum_i { \sum_{\nu} { c_{i\mu}c_{i\nu} S_{\mu\nu} \delta(\varepsilon - \varepsilon_i) } }\]

which is the projected density of states (pDOS) over orbital \(\mu\). Note that when the basis is not orthogonal, the sum over \(\mu\) and the overlap factor are needed.

The tools available in Siesta to compute and process the DOS and pDOS are discussed in this how-to.

COOP and COHP curves

Yes!. We can take out another sum in the above equation:

\[g_{\mu\nu} (\varepsilon) = \sum_i { c_{i\mu}c_{i\nu} S_{\mu\nu} \delta(\varepsilon - \varepsilon_i) }\]

This, if we follow the chain of reasoning, must be the contribution of the pair of orbitals \(\mu, \nu\), to the electron charge “at some energy”. If we had started from the formal decomposition of the band-structure energy, we would have arrived at:

\[h_{\mu\nu} (\varepsilon) = \sum_i { c_{i\mu}c_{i\nu} H_{\mu\nu} \delta(\varepsilon - \varepsilon_i) }\]

which similarly would be the pair contribution to the (band structure) energy.

The above two equations define what are known, respectively, as the Crystal Orbital Overlap Population (COOP) and Crystal Orbital Hamilton Population (COHP) curves. (Curves because they are functions of energy).

It turns out that one can obtain a lot of insight into the chemical bonding from these curves.

For background and examples, see the LOBSTER site for more information. LOBSTER is a program created by the group of Richard Dronskowski, one of the pioneers of this kind of bonding analysis in solid-state chemistry.

In this tutorial we offer a couple of examples to show the implementation in Siesta of this analysis machinery:

Fat-bands: orbital-projected band-structure

This technique is not (directly) related to the chain of equations discussed above, but is quite useful to get a visual idea of the character of the bands.

Analysis of the spin texture in reciprocal space

This is more advanced but provides special insight.