Methodology

The pCOHP and the pCOBI

The projected crystal orbital Hamilton population (pCOHP) and the projected crystal orbital bond index (pCOBI) are both local descriptors of chemical bonding defined by the projection of the Kohn-Sham eigenstates (as calculated via a prior DFT calculation) onto a chosen set of localised basis functions [1][2]

\[\begin{split}\begin{align} \mathrm{pCOHP}_{\alpha\beta}(E) &= H_{\alpha\beta}\sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\beta}_{nk}\right]\cdot\delta(E - \epsilon_{nk}) \\ \mathrm{pCOBI}_{\alpha\beta}(E) &= P_{\alpha\beta}\sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\beta}_{nk}\right]\cdot\delta(E - \epsilon_{nk}). \end{align}\end{split}\]

In both of the equations above, each \(C^{\alpha}_{nk}\) is the coefficient resulting from the projection of the Kohn-Sham eigenstate \(\ket{\psi_{nk}}\) onto the localised basis function \(\ket{\phi_{\alpha}}\). Together with \(\delta(E - \epsilon_{nk})\), which is the density of states arising from \(\ket{\psi_{nk}}\), everything within the summation in both equations is referred to as the DOS matrix arising from local basis states \(\ket{\phi_{\alpha}}\) and \(\ket{\phi_{\beta}}\):

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

In the pCOHP, the sum over this DOS matrix is weighted by the corresponding element of the local basis Hamiltonian \(H_{\alpha\beta}\), whilst in the pCOBI it is weighted by \(P_{\alpha\beta}\), the relevant element of the local basis density matrix.

Loosely speaking, the pCOHP is thought to be correlated with bond strength, whilst the pCOBI is thought to be correlated with bond order (or rather, the contribution to these two metrics by a given pair of basis functions). This is most clearly expressed by their integrals up to the Fermi level,

\[\begin{split}\begin{align} \mathrm{ICOHP}_{\alpha\beta} &= \int^{E_{\mathrm{F}}}_{-\infty} \mathrm{d}E\,\mathrm{COHP}_{\alpha\beta}(E) \\ \mathrm{ICOBI}_{\alpha\beta} &= \int^{E_{\mathrm{F}}}_{-\infty} \mathrm{d}E\,\mathrm{COBI}_{\alpha\beta}(E), \end{align}\end{split}\]

which have often been used as a quantitative measure of these two characteristics.

The spilling factor

Perhaps the most serious problem with calculating pCOHPs and pCOBIs in the manner detailed above is that the one requires a suitable set of localised basis functions \(\ket{\phi_{\alpha}}\) that, when appropriately combined, make for a sufficiently accurate representation of the original Kohn-Sham eigenstates. If the local basis is not sufficiently representative of the canonical Bloch states, then there will be a net loss of information after the projection and the resulting pCOHPs and pCOBIs will be correspondingly less accurate. This potential loss of information is quantified by the spilling factor [3]

\[S = \frac{1}{N_{k}}\frac{1}{N_{b}}\sum_{nk} 1 - \sum_{\alpha}|\braket{\psi_{nk}|\phi_{\alpha}}|^{2},\]

which takes values between 0 and 1. If the local basis spans the same Hilbert space as the Kohn-Sham states, then \(S = 0\), whilst \(S = 1\) indicates that the two bases are orthogonal to one another.

Wannier functions

Very briefly, Wannier functions are a localised basis obtained from a set of Kohn-Sham eigenstates via a unitary transformation [4]. For a set of \(J\) energetically isolated bands (i.e., a manifold of bands separated from all other bands by an energy gap everywhere in the Brillouin zone), this can be written as

\[\begin{split}\begin{align} \ket{w_{\alpha}} &= \sum_{k}\exp[-ik\cdot R]\sum^{J}_{m} U^{k}_{m\alpha}\ket{\psi_{mk}} \\ \ket{w_{\alpha}} &= \ket{w_{iR}}, \end{align}\end{split}\]

where \(i\) is a band-like index (note that within the exponential term, \(i = \sqrt{-1}\)), \(R\) is a Bravais lattice vector specifying the unit cell in which \(\ket{w_{iR}}\) resides (\(\alpha\) combines these two indices) and each \(U^{k}\) is a unitary matrix. Owing to the fact that the Wannier functions are obtained via unitary transformations, the trace over any one-particle operator is the same whether one chooses to use the Bloch basis or the Wannier basis:

\[\sum_{k}\ket{\psi_{nk}}\bra{\psi_{nk}} = \sum_{R}\ket{w_{nR}}\bra{w_{nR}}.\]

The spilling factor can therefore be written as

\[S = \frac{1}{N_{k}}\frac{1}{N_{b}}\sum_{nk} 1 - \sum_{m}|\braket{\psi_{nk}|\psi_{mk}}|^{2},\]

which due to the orthonormality of the Kohn-Sham states is simply:

\[S = \frac{1}{N_{k}}\frac{1}{N_{b}}\sum_{nk} 1 - \sum_{m}\delta_{mn} = 0.\]

In the context of chemical bonding descriptors, Wannier functions therefore represent an ideal localised basis for the computation of the pCOHP and pCOBI or alternatively the WOHP [5] and WOBI (W for Wannier):

\[\begin{split}\begin{align} \mathrm{WOHP}^{R}_{\alpha\beta}(E) &= -H_{\alpha\beta}^{R}\sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\beta}_{nk}\right]\cdot\delta(E - \epsilon_{nk}) \\ \mathrm{WOBI}^{R}_{\alpha\beta}(E) &= P_{\alpha\beta}^{R}\sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\beta}_{nk}\right]\cdot\delta(E - \epsilon_{nk}). \end{align}\end{split}\]

Note

We define the WOHP with the opposite sign to the pCOHP so as to ensure that a positive WOHP indicates bonding and a negative WOHP indicates antibonding (thus matching the WOBI in this respect).

where \(R = R_{2} - R_{1}\) is an extra index accounting for the fact that one could techincally compute the WOHP or WOBI between Wannier functions located in a variety of different unit cells. The coefficients used to build the DOS matrix are easily obtained from the unitary matrices used to define the Wannier functions:

\[C^{\alpha}_{nk} = \exp[ik \cdot R]\left(U^{k}_{n\alpha}\right)^{*}.\]

The pCOOP and the WOOP/pDOS

As well as the pCOHP and the pCOBI, there is also a third descriptor that has traditionally been used to assess the bonding/antibonding character of various interactions: the projected crystal orbital overlap population or pCOOP:

\[\mathrm{pCOOP}_{\alpha\beta}(E) = S_{\alpha\beta}D_{\alpha\beta}(E) = S_{\alpha\beta}\sum_{nk}\mathrm{Re}\left[\left(C^{\alpha}_{nk}\right)^{*}C^{\beta}_{nk}\right]\cdot\delta(E - \epsilon_{nk}),\]

where \(S\) is the overlap matrix. If we use Wannier functions as our local basis, then \(S = I\) and the WOOP (Wannier orbital overlap population) is non-zero only for on-site interactions [5]:

\[\mathrm{WOOP}_{\alpha\alpha}(E) = D_{\alpha\alpha}(E) = \mathrm{pDOS}_{\alpha}(E),\]

in which case we can also think of it as \(\mathrm{pDOS}_{\alpha}\): the density of states projected onto a specific Wannier function \(\ket{w_{\alpha}}\).

Caveats

Whilst Wannier functions derived from energetically isolated bands are guaranteed to span the same Hilbert space as the Kohn-Sham eigenstates, this does not mean that individual WOHPs and WOBIs are uniquely defined. In general, Wannier functions are strongly non-unique: so long as each \(U^{k}\) remains unitary, the spilling factor is strictly 0, but the resulting Wannier functions may have wildly different centres, spreads and shapes. In addition, Wannier bases derived from groups of entangled bands (those that are not separated by an energy gap from all other bands everywhere in the Brillouin zone) are no longer guaranteed to span exactly the same space as the original Kohn-Sham states, thus potentially suffering from the same problems as pre-defined basis sets of atomic or pseudo-atomic orbitals.

There exists well-estalished methods for mitigating both of the concerns raised above. The non-uniqeuness of Wannier functions can be circumvented by minimising their spread with respect to the unitary matrices \(U^{k}\), thus producing so-called “Maximally-Localised Wannier Functions” or MLWFs [6][7]. Such Wannier functions tend to be atom-centred and to take on shapes that in many cases match our chemical intution for the system at hand. Maximal localisation can also be applied to a manifold of entangled bands, although a separate “disentanglement” step is required in this case, unless one chooses to utilise the SCDM method [8][9].

In general, obtaining “good” Wannier functions for non-trivial systems is not always easy, although developments in recent years are slowly making this a more reliable and automatable process [10][11]. As a result, the two most obvious use cases for pengwann are as follows:

1. You have already obtained MLWFs for your system (for some other purpose) and would like to characterise the local bonding.

2. You have attempted to characterise the bonding via another tool such as LOBSTER, but you find that the results are not satisfactory.

References

---

[1]

Volker L. Deringer, Andrei L. Tchougréeff, and Richard Dronskowski. Crystal orbital Hamilton population (COHP) analysis as projected from plane-wave basis sets. J. Phys. Chem. A, 115(21):5461–5466, 2011. PMID: 21548594. URL: https://doi.org/10.1021/jp202489s, arXiv:https://doi.org/10.1021/jp202489s, doi:10.1021/jp202489s.

[2]

Peter C. Müller, Christina Ertural, Jan Hempelmann, and Richard Dronskowski. Crystal orbital bond index: covalent bond orders in solids. J. Phys. Chem. C, 125(14):7959–7970, 2021. URL: https://doi.org/10.1021/acs.jpcc.1c00718, arXiv:https://doi.org/10.1021/acs.jpcc.1c00718, doi:10.1021/acs.jpcc.1c00718.

[3]

Daniel Sanchez-Portal, Emilio Artacho, and Jose M Soler. Projection of plane-wave calculations into atomic orbitals. Solid State Commun., 95(10):685–690, 1995. URL: https://www.sciencedirect.com/science/article/pii/003810989500341X, doi:10.1016/0038-1098(95)00341-X.

[4]

Gregory H. Wannier. The structure of electronic excitation levels in insulating crystals. Phys. Rev., 52:191–197, Aug 1937. URL: https://link.aps.org/doi/10.1103/PhysRev.52.191, doi:10.1103/PhysRev.52.191.

[5] (1,2)

Sudipta Kundu, Satadeep Bhattacharjee, Seung-Cheol Lee, and Manish Jain. Population analysis with Wannier orbitals. J. Chem. Phys., 154(10):104111, 03 2021. URL: https://doi.org/10.1063/5.0032605, arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0032605/13469861/104111\_1\_online.pdf, doi:10.1063/5.0032605.

[6]

Nicola Marzari and David Vanderbilt. Maximally localized generalized Wannier functions for composite energy bands. Phys. Rev. B, 56:12847–12865, Nov 1997. URL: https://link.aps.org/doi/10.1103/PhysRevB.56.12847, doi:10.1103/PhysRevB.56.12847.

[7]

Ivo Souza, Nicola Marzari, and David Vanderbilt. Maximally localized Wannier functions for entangled energy bands. Phys. Rev. B, 65:035109, Dec 2001. URL: https://link.aps.org/doi/10.1103/PhysRevB.65.035109, doi:10.1103/PhysRevB.65.035109.

[8]

Anil Damle, Lin Lin, and Lexing Ying. Compressed representation of Kohn–Sham orbitals via selected columns of the density matrix. J. Chem. Theory Comput., 11(4):1463–1469, 2015. PMID: 26574357. URL: https://doi.org/10.1021/ct500985f, arXiv:https://doi.org/10.1021/ct500985f, doi:10.1021/ct500985f.

[9]

Anil Damle, Lin Lin, and Lexing Ying. SCDM-k: localized orbitals for solids via selected columns of the density matrix. J. Comp. Phys., 334:1–15, 2017. URL: https://www.sciencedirect.com/science/article/pii/S0021999116307215, doi:https://doi.org/10.1016/j.jcp.2016.12.053.

[10]

Valerio Vitale, Giovanni Pizzi, Antimo Marrazzo, Jonathan R. Yates, Nicola Marzari, and Arash A. Mostofi. Automated high-throughput Wannierisation. npj Comput. Mater., 6(1):66, Jun 2020. URL: https://doi.org/10.1038/s41524-020-0312-y, doi:10.1038/s41524-020-0312-y.

[11]

Junfeng Qiao, Giovanni Pizzi, and Nicola Marzari. Projectability disentanglement for accurate and automated electronic-structure Hamiltonians. npj Comput. Mater., 9(1):208, Nov 2023. URL: https://doi.org/10.1038/s41524-023-01146-w, doi:10.1038/s41524-023-01146-w.