LOBSTER: A tool to extract chemical bonding from plane-wave based DFT.

The computer program LOBSTER (Local Orbital Basis Suite Towards Electronic-Structure Reconstruction) enables chemical-bonding analysis based on periodic plane-wave (PAW) density-functional theory (DFT) output and is applicable to a wide range of first-principles simulations in solid-state and materials chemistry. LOBSTER incorporates analytic projection routines described previously in this very journal [J. Comput. Chem. 2013, 34, 2557] and offers improved functionality. It calculates, among others, atom-projected densities of states (pDOS), projected crystal orbital Hamilton population (pCOHP) curves, and the recently introduced bond-weighted distribution function (BWDF). The software is offered free-of-charge for non-commercial research.

Introduction

Methods for electron partitioning in molecules have been around in quantum chemistry since Mulliken's ingenious approach for assigning electrons to atoms and bonds.1 These models and concepts are likewise helpful for periodic systems, so an analogous scheme was introduced within non‐variational extended Hückel (eH) theory2 and dubbed Crystal Orbital Overlap Population (COOP); pioneered in the 1980s, it proved powerful ever since.3 With the advent of variational density‐functional theory (DFT), the Crystal Orbital Hamilton Population (COHP) scheme was suggested, which partitions energies rather than electrons but otherwise resembles COOP in that it allows to extract chemical interactions between atoms from band‐structure calculations.4 COHPs have been implemented and widely used within TB‐LMTO‐ASA theory,5 which is DFT‐type but shares with eH the use of localized basis sets for periodic solids.

As of today, many condensed‐matter quantum‐mechanical codes employ plane waves (PW), which naturally (and effectively) fulfill Bloch's theorem but are delocalized by their very nature, making bond‐analytical approaches such as COOP and COHP unavailable in PW frameworks. Nonetheless, there are ways to transfer PW to localized functions using projection schemes as pioneered by Sánchez‐Portal et al.6 For large atomic numbers, however, PW become impractical for the near‐core regions, so the success of the pseudopotential (PP) approach7 in computational materials science is easily understandable. Nowadays, Blöchl's projector‐augmented wave (PAW) method is the most powerful of the PP descendants.8 To project PAW functions8 onto localized orbitals (say, of the Slater type), we have recently developed an analytical formalism9 to apply bond‐analytic tools even though the system was brought to self‐consistency in a PW basis; other bond‐analytical approaches10 exist as well. Our technique makes COOP or COHP analyses feasible beyond densely packed systems9 (such as intermetallics11) and can, other than before, be seamlessly applied to scenarios such as molecular crystals12 or even amorphous matter.13

To facilitate chemical‐bonding analyses and other methods for a multitude of systems and applications, we are offering the LOBSTER (Local Orbital Basis Suite Towards Electronic‐Structure Reconstruction) software free‐of‐charge for non‐commercial purposes at http://www.cohp.de. In this article, we describe recently developed methodology and functionality added to LOBSTER. A number of illustrative applications are presented, and directions for further reading are given.

Methods

Before describing recent developments which have found their way into LOBSTER, we refer the reader to our initial publication9 for a more comprehensive account of the underlying formalisms, and to the manual shipped with the code for any practical questions. A (simplified) scheme of what LOBSTER does is presented in Figure 1.

Figure 1

Overview of LOBSTER's functional principle: a quantum‐chemical system, characterized by its one‐electron (Bloch) wavefunctions and the according eigenvalues (band energies), has been brought to self‐consistency using some plane‐wave DFT program. A local auxiliary basis is then selected to determine the overlap matrix and the transfer matrix between the delocalized and localized representations. From those, the projected coefficient and Hamiltonian matrices and , respectively, are accessible, which allow for various bond‐analytic tools. The LOBSTER logo is copyrighted by the Chair of Solid‐State and Quantum Chemistry at RWTH Aachen University.

Local basis sets

The initial step in projection is finding a suitable choice of local auxiliary basis functions. For reasons of simple chemical interpretation, LOBSTER employs minimal basis sets that nonetheless carry the correct nodal behavior in the core region, which is necessary to fit PAW wavefunctions. LOBSTER first came with contracted primitive Slater‐type orbitals (STOs) fitted to atomic functions,14 a reasonable choice for post‐processing bonding information. There are also systems, however, where the bonding situation requires additional basis functions which deviate from those of the corresponding free atoms. In elemental beryllium, for example, its 2p levels are unoccupied in the free atom. For bulk Be, however, the 2p levels do get involved into bonding and define the metallic character, so the Be basis set must also include 2p functions.

For demonstration, let's look at the high‐temperature phase,15 body‐centered cubic beryllium, ‐Be. Its electronic structure was calculated with ABINIT employing the JTH atomic datasets16 and the GGA‐PBE parametrization for exchange and correlation.17 On the LOBSTER side, the original basis set[14a] and its basis functions (1s and 2s) somehow manage to reconstruct the PAW electronic structure but with an unacceptable absolute charge spilling of roughly 19% (see below for definitions). For analysis, the differences between the original and projected wavefunctions were calculated and an isosurface at 65% of the maximum resulting density was plotted for several bands at . The fourth band showed enormous deviations (Fig. 2, left) because the basis lacks an orbital of p‐symmetry, as reasoned before. Adding a 2p function strongly reduces the absolute charge spilling to 1.73%. For comparison, the 65% density‐difference maxima decrease by two orders of magnitude (Fig. 2, right). While the functions were added by fitting VASP data, they turn out to be general enough to easily fit other PAW wavefunctions, for example, those calculated by ABINIT, too.

Figure 2

Isosurfaces (in Å−3) at 65% of the differences between the ABINIT‐based PAW densities and the LOBSTER‐projected densities for the fourth band of ‐Be at . On the left side, the basis contains only 1s and 2s functions as given by Bunge et al., whereas on the right this basis was enlarged by 2p functions. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

Hence, free‐atom calculations in large supercells have been performed for all elements up to (curium) using GGA‐PBE17 as implemented in VASP.18 In nearly all cases up to (mercury) did the new basis functions match the previously given ones well, and they allowed us to numerically fit and add missing (polarization) functions. Obviously, these new functions had to be orthogonalized with regard to the existing functions of the same l azimuthal quantum number to enlarge the basis sets already available in LOBSTER. In the next step, visual evaluation of the calculated PAW atomic orbitals yielded wavefunctions of the desired symmetry and shape; hence they were corroborated as a valid basis choice.

While the basis functions are aligned with the Cartesian axes by default, LOBSTER 2.0.0 supports user‐defined rotations of the basis functions as has been described recently;19 this new feature can be especially useful when isolated, orbital‐wise interactions must be studied.

Improved measures for projection quality

Before analyzing the projected wavefunctions, one must ensure that the auxiliary basis suffices. To do so, Sánchez‐Portal et al. introduced the so‐called “spilling” and “charge‐spilling” criteria,6 b which measure the percentage of electronic density lost during projection; this approach was also used in former LOBSTER versions.9

Nonetheless, in a localized basis, unwanted effects such as the basis‐set superposition error20 or counterintuitive orbital mixing21 can lead to projected wavefunctions with a norm artificially larger than unity; note that the original spilling criterion correctly assumes that the norm of a projected wavefunction is bound to unity. If this condition is broken, averaging the spilling over multiple bands and points may lead to error cancellation such that the projection quality looks better than it actually is. To counteract, LOBSTER 2.0.0 comes with an improved definition which we dub “absolute spilling” and “absolute charge spilling”, , defined analogously to its predecessor,6 b but averaging absolute values: with where and are the coefficient and overlap matrices, respectively, and denotes the normalized point weights. The absolute charge spilling, , is collected only over occupied bands, that is, those with nonzero occupation numbers, .

Another way to assess the deviation is given by the root‐mean‐square error (RMS) of projected wavefunctions. If the PW part of the PAW functions is given on a reciprocal grid at values of , one may define the RMS of the projection (RMSp) by comparing the projected LCAO function to its PAW reference :

The sum runs over all vectors at each band j in the plane‐wave basis. Using the reciprocal representation of the wavefunctions enables us to rewrite the difference in the former equation:

All of these expressions are either known directly from the PAW calculation (like the plane‐wave coefficients and so‐called wavefunction characters , evaluated during the projection by LOBSTER anyway [such as the LCAO coefficients and the Fourier‐transforms of the local basis functions ], or they can be obtained by a Fourier–Bessel transform where is a shorthand notation for the difference between the all‐electron and pseudo‐space partial‐waves in the PAW method and . designates the spherical Bessel function, and is a spherical harmonic.

In contrast to the originally defined spilling, the RMSp method is bound to zero independent of assumptions, which makes it a well‐suited optimization criterion. If desired, one may normalize RMSp to the range of reference data, viz. , for example, when comparing results from primitive unit cells to those from supercell models, a rather practical real‐world scenario.

Orthonormalization

As stated before, the targeted analytic methods are bound to minimal basis sets on purpose and hence prevent using sophisticated multi‐ basis sets easily. Consequently, we apply Löwdin's symmetric orthonormalization (LSO) to the projected wavefunctions.9 Even if significant (but comparable) amounts of charge are lost around every atom, this technique was found to ensure properly projected densities of states (pDOS), a crucial ingredient for other bond‐analytical tools in the sequel.

Within LOBSTER 2.0.0, LSO is also applied to the basis functions themselves. We note that traditional COHP analysis by means of TB‐LMTO‐ASA theory5 works with an intrinsically orthogonal basis set which overlaps due to the atomic‐spheres‐approximation: likewise, the basis functions in LOBSTER do overlap. To improve correspondence between traditional COHP and its projected analogue, the projected Hamiltonian matrix is now reconstructed after likewise applying LSO to the basis functions, yielding where designates the coefficient matrix within the orthogonalized set of basis functions. In contrast, all other bond‐analytic tools such as projected COOP still use since the overlap populations would be rendered meaningless in an orthogonal basis set.

Visualization

Examining the causes of an imperfect projection is a non‐trivial task but may be significantly simplified by visual inspection. To do so, the internal development version of LOBSTER writes the values of both projected and PAW wavefunctions on a user‐defined, linearly equidistant grid which can either be an arbitrarily oriented line or a three‐dimensional grid within a cuboid, both bounded by the unit cell. This enables density‐difference plots as shown in Figure 2, but can also be used to examine the signed wavefunctions directly.

Technical aspects

LOBSTER is written in modern, object‐oriented C++ and uses the famous Boost library22 for various algorithms and concepts for object or memory management, as well as their organization. Even though many of them were incorporated into the C++11 standard, Boost is still a valuable asset for mathematical special functions and interaction with the operating system. Furthermore, LOBSTER has evolved into a multiplatform tool, supporting Linux, Windows and OS X. Wherever possible, it employs matrix or vector algebra to employ the data structures and algorithms, for example, for matrix decompositions, provided by the highly efficient Eigen library.23 Results of computationally expensive but repeatedly evaluated functions are cached internally. LOBSTER is parallelized using OpenMP, still making efficient use of its internal caches through multiple‐readers/single‐writer lock patterns. To optimize CPU usage, LOBSTER uses memory mapped I/O to read large chunks of input data from the file system when available and beneficial.

At present, LOBSTER 2.0.0 processes and analyzes PAW results from two third‐party codes, VASP18 and ABINIT,24 but further interfaces are possible. Due to lack of direct tabulation in the case of ABINIT, the projector functions must be transformed to reciprocal space once using the Fourier–Bessel transform, as stated above.

Program Features

Once the coefficient matrix has been reconstructed (Fig. 1), LOBSTER can readily calculate pDOS and pCOOP by space integration.9 Additionally, it writes their energy‐integrated counterpart IpDOS (which yields the total number of electrons of the respective atoms, Mulliken's gross population) and IpCOOP (Mulliken's overlap population). Based on these IpCOOP values, bonding analysis can be applied even to amorphous structures by means of the recently proposed bond‐weighted distribution function (BWDF).13 The coefficient matrices have also been used directly to analyze orbital mixing (or hybridization in the physicists' language).25 Reconstructing the Hamiltonian matrix in a second step facilitates pCOHP analysis.9, 26 Energy integration up to the Fermi level yields IpCOHP, which might likewise serve as a bond‐weighting indicator for BWDF.13 For further detail on each specific method, we redirect to the original literature.

Applications

Since its initial publication, LOBSTER has found diverse applications. Being interfaced to state‐of‐the‐art DFT codes such as VASP, it allows to process output from modern methods such as hybrid‐DFT results.27 LOBSTER has begun to play its part in surface chemistry: exploring oxide catalysts,28 square‐planar carbon at transition‐metal surfaces,29 or local structural fragments at quartz‐type GeO2 surfaces.30 Less‐than‐densely packed three‐dimensional (3D) networks have been of interest as well: complex clathrate structures,31 hydrogen bonding in molecular crystals,12 or the stability ranking of metal azide polymorphs.32

We round out this article by presenting two representative applications from fields of current research interest.

Atomic motion in phase‐change materials

Phase‐change materials (PCMs) can be switched between crystalline and amorphous phases, and thus be used to encode “ones” and “zeroes” in digital data storage.33 Atomistic simulations, generally a cornerstone of PCM research,34 have recently been concerned with transition pathways and atomic motion in crystalline PCMs.35

The prototypical PCM germanium telluride (GeTe) has been studied using COHPs more than a decade ago using TB‐LMTO‐ASA theory.36 With the new projection techniques at hand, we may now explore not only its crystalline and amorphous phases,13 but also the formation and diffusion of vacancies on the crystalline (distorted rocksalt‐type) lattice, as visualized in Figure 3 (top). The well‐known presence of antibonding interactions in PCMs37 is also seen in LOBSTER output (Fig. 3, bottom). While these regions (–pCOHP < 0) are there both for the moving atom (red) and away from the transition state (gray), they critically depend on the local environment: the antibonding peaks are much more dominant in the transition‐state geometry (Fig. 3, bottom left). This reflects the moving Ge atom which has to “squeeze” through an octahedral face to leave its ground‐state position, which costs over 100 kJ/mol according to DFT simulations.35 a Further away, on the contrary, the bonding situation seems relatively unperturbed (Fig. 3, bottom right) and compares well to what is found in the ground‐state structure. The bonding analysis also allows us to rationalize previously made observations regarding Fermi level tuning: lowering it makes the transition state (‡) less unfavorably bonded, and thus reduces the activation energy for the Ge hop significantly.35 a, 38

Figure 3

LOBSTER analysis of a diffusion pathway through crystalline GeTe, as originally mapped out using nudged‐elastic‐band (NEB) theory35 a and previously analyzed by pCOOP analysis in the thesis of one of us.38 Top: structural drawing of the supercell setup, with only selected atoms shown for clarity.38 A germanium atom jumps from one octahedron into an adjacent one; a second one further away serves as reference. Bottom: pCOHP analysis for the sum of the three short Ge–Te bonds shown, respectively. Energy is shifted so that the Fermi level equals zero.

The orbital origins of ferromagnetism, revisited: from ‐iron to Ru monolayers

Ferromagnetism of ‐Fe can be chemically understood as a consequence of orbital interactions which has previously been rationalized at the hand of COHP analysis,39 and this concept later served as a predictive guideline for the design of more complex magnetic materials.11, 40 Figure 4 (left) presents a COHP analysis based on a traditional TB‐LMTO‐ASA calculation at the spin‐restricted (non‐magnetic, NM) GGA level showing occupied antibonding levels at the Fermi level which destabilize the system. The middle and right of Figure 4 offer analogous analyses, but employ the methods and frameworks given in this work. In the non‐magnetic case, the VASP/LOBSTER combination recovers what has been known before. By switching on spin‐polarization, the majority (α) electrons lower in energy (and its associated spin orbitals spatially contract) while the minority (β) electrons increase energetically (and its orbitals expand);39 as a result, antibonding states are diminished and the chemical bonding strengthens.

Figure 4

Bonding analysis of ‐iron using (left) COHP as implemented in TB‐LMTO‐ASA in a non‐magnetic (NM) setup and (middle) using pCOHP based on PAW results by VASP processed with LOBSTER. By allowing for spin polarization (SP, right), the resulting exchange splitting affects the chemical bonding between the Fe atoms which becomes stronger. COHP and pCOHPs are given as the sum over all symmetry‐equivalent bonds in the unit cell. Energy is shifted so that the Fermi level equals zero. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

While the above is merely validation, the new method (in contrast to TB‐LMTO‐ASA theory) can easily handle “open” systems such as two‐dimensionally extended surface structures. Magnetism in such systems has been explored earlier, and one of them is shown in Figure 5 (left): a monolayer of ruthenium atoms supported on a slab of Ag. As originally predicted by Blügel,41 ruthenium becomes ferromagnetic in this configuration. Figure 5 (right) now delves into the chemical‐bonding nature again, and it suggests an explanation that is principally analogous to the ‐iron case: namely significant exchange splitting and strengthening of the Ru–Ru bonds while becoming ferromagnetic, now so easily rationalized using LOBSTER.

Figure 5

Spin polarization and chemical bonding in a two‐dimensional Ru sheet supported on an Ag(001) surface. The left‐hand side shows the supercell setup, and a large “vacuum” area is clearly visible, as commonly used in PW based DFT simulations of surfaces and nanomaterials. The right‐hand side shows the LOBSTER‐computed pCOHP curve for a single nearest‐neighbor Ru–Ru contact in the spin‐polarized case; compare with Figure 4 (right). Energy is shifted so that the Fermi level equals zero. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

Conclusions

We have presented new developments in the LOBSTER software for chemical‐bonding analysis. LOBSTER processes delocalized PAW wavefunctions calculated with VASP or ABINIT and performs projection into an auxiliary LCAO basis, which makes bond‐analytic tools such as pDOS, pCOOP, and pCOHP accessible for state‐of‐the‐art plane‐wave based PAW simulations.

To reliably assess the quality of the projection, we here introduced two modified criteria, dubbed absolute (charge) spilling and root‐mean‐square of the projection (RMSp). Additionally, visual evaluation of either the PAW or the projected wavefunctions has been demonstrated. A new and improved basis set available in LOBSTER has also been described. Finally, to improve correspondence to traditional COHP analysis based on LMTO theory, we now also apply Löwdin's symmetric orthogonalization to the basis functions.