Automated Bonding Analysis with Crystal Orbital Hamilton Populations
Understanding crystalline structures based on their chemical bonding is growing in importance. In this context, chemical bonding can be studied with the Crystal Orbital Hamilton Population (COHP), allowing to quantify interatomic bond strength. Here we present a new set of tools to automate the calculation of COHP and analyze the results. We use the program packages VASP and LOBSTER and the Python packages atomate and pymatgen. The analysis produced by our tools includes plots, a textual description, and key data in machine-readable format. To illustrate those capabilities, we have selected simple test compounds (NaCl, GaN), the oxynitrides BaTaO2N, CaTaO2N, and SrTaO2N, and the thermoelectric material Yb14Mn1Sb11. We show correlations between bond strengths and stabilities in the oxynitrides, as well as the influence of the Mn-Sb bonds on the magnetism in Yb14Mn1Sb11. Our contribution enables high-throughput bonding analysis and will facilitate the use of bonding information for machine learning studies.
Introduction<!>The workflow<!>Automatic Test of Best Basis and Projections<!>Two illustrative examples: NaCl and GaN<!>Two advanced examples: oxynitrides and Yb14Mg1Sb11<!>Conclusion<!>Methods
<p>The concept of chemical bonding plays a major role in understanding both molecules and crystals and their properties. [1,2] Chemical bonding in crystals can be studied using various techniques based on electronic structure theory methods. These techniques are usually divided into density-based and orbital-based methods. For example, a widely used densitybased method is Bader analysis, which is based on electron density topology analysis. [3] Bader analysis has also been used in a high-throughput manner. [4][5][6] In this paper, however, we will focus on a more quantum-chemical orbital-based bonding analysis, which has been shown to be successful in understanding crystal structures in the past, [2,[7][8][9] and its automation. These methods have the advantage that bonding and antibonding contributions directly fall off from the orbital phases such that one arrives at a picture that better relates to the chemist's orbital-based understanding of atoms and molecules. [8] An important technique in this context is Crystal Orbital Hamilton Population (COHP), [10] a bond-weighted density of states that provides information about bond energy and covalency in crystals. [7] There are other bond indicators not based on the Hamilton matrix but on the overlap matrix or bond order (Crystal Orbital Overlap Populations, COOP and Crystal Orbital Bond Index, COBI). [11] The first one relates to bond strengths, the second one to bond orders. In addition, there are also orbitalbased indicators for the ionicity of compounds. [11,12] COHP has been used for about three decades to study all kinds of compounds, for example, phase change materials, [13] intermetallics, [14] magnetic materials, [15] and many other classes of materials. Nowadays, these populations can also be derived according to electronic structure theory calculations using the projector augmented wave method. This is done by projecting the original delocalized basis based on plane waves onto atomic orbitals (e.g., Slater orbitals), as implemented, for example, in the computer program LOBSTER. [16][17][18][19] LOBSTER is also capable of computing COHPs, COOPs, COBIs, and many other similar properties based on the information resulting from the projection. In addition, Mulliken and Löwdin charges and corresponding Madelung energies can be calculated. [11,20] This procedure is illustrated in Figure 1. Here we use VASP as the code for our electronic structure theory calculations.</p><p>To date, this type of projection has not been implemented in any of the standard codes for electronic structure theory based on plane waves (e.g., VASP, [21][22][23] Abinit, [24] or Quantum Espresso [25] ). Therefore, this type of bonding analysis is a multi-step process involving at least two programs. New users not only have to use and learn non-standard settings for the electronic structure theory codes, they must also learn to use LOBSTER. Naturally, the users should also be aware of fundamental quantum-chemical concepts not necessarily taught in courses of quantum mechanics or computational materials science. Hence, this usually involves a lot of training for first-time users of the programs, a lot of manual work, possible errors and problems in interpretation. This limits the use of this bonding analysis tool for high-throughput calculations, although there have been some studies where LOBSTER has been used for a larger number of compounds. [26] Workflow-managing codes nowadays provide ways to automate such complex calculations. [6,27,28] This has been done in the past for several types of standard calculations (e.g., magnetic groundstate calculations, defect calculations, [6,[29][30][31][32] ). This automation is also directly related to the need for high-throughput calculations and accurate computational data in materials science for machine learning and data-driven material discovery. [33][34][35] Figure 1. Illustration of the process of bonding analysis using the VASP and LOBSTER tools. VASP is used to calculate a wave function using the projector-augmented wave method. This calculation depends on the basis function chosen in the later LOBSTER run. LOBSTER then projects this wave function onto a basis set of atomic orbitals that must be adapted to the system at hand, corresponding to a unitary transformation between a totally delocalized and a localized representation. Then, LOBSTER calculates various populations and charges based on these populations, which allow to assess the bonding situation in the compound. This is a multi-step process whose steps are interdependent. We have now developed automated tools to perform all the steps of this procedure and to ensure the compatibility of the different steps.</p><p>For example, high-throughput bonding analysis could provide further insight into important interactions within crystal structures and provide opportunities to determine coordination environments in crystal structures not only based on geometric constraints but also on electronic structure. [36,37] This has the advantage that the relevance of neighbors can be evaluated based on the integrated COHP as a measure of a bond energy and hence a covalent bond strength. In purely geometric approaches, the strength of covalent bonds can only be evaluated indirectly using limits on bond lengths. Moreover, some of these bond properties could serve as descriptors in machine learning studies since the underlying bonding situation is expected to be causative for many materials properties. Recently, two such studies have been published. [26,38] Here, we present an implementation of bonding analysis using the programs VASP and LOBSTER in a Python-based workflow. We have automated this process using the programs VASP, LOBSTER, [7,[16][17][18][19] pymatgen, [39] custodian, [39] and atomate. [28] This means that the entire bonding analysis workflow, leading to an analysis text and a summary of key properties, can be started with very little Python code. In addition to the computational workflow, we have also automated the analysis of LOBSTER results. Without this automation, users must do this by hand, and there is currently no unified strategy for doing so, unless the user is properly trained. Therefore, automating this process will also help in standardizing the analysis of LOBSTER results. This automation of the whole process will hopefully lead to fewer errors in the future and less training being required. We also expect that the automated analysis will contribute to a much wider use of the tools and serve as a starting point for more detailed and individual bonding analysis by the user. In addition, high-throughput studies should now be in reach for non-expert LOBSTER users.</p><p>To illustrate the capabilities of the code, we have tested it on a range of binary and more complex materials: we use the binaries NaCl, GaN, then the oxynitrides CaTaO2N, BaTaO2N, SrTaO2N, and finally Yb14Mn1Sb11, where the bonding analysis is performed fully automatically. In the case of the oxynitrides, we also show how easily correlations with other material properties can be derived based on the developed tools. We test the correlation between the total energy of the systems and the strongest integrated COHP upon varying the anion order in these perovskite systems. In addition, Yb14Mn1Sb11 illustrates the potential of the code to identify important interactions for magnetic properties and potential bonding-based driving forces for properties in materials.</p><!><p>The implemented workflow is shown in Figure 2.</p><p>Figure 2. This graphic illustrates the progress that has been made in this work. Instead of creating all input files through custom scripts or by hand and starting and evaluating calculations by hand, there is now a fully automated workflow for the entire bonding analysis process using the VASP and LOBSTER programs. This workflow can automatically create all input files, start and monitor calculations. We have also developed tools to evaluate all the results and create automatic output texts and dictionaries of the most important bonding properties, which can be used to populate databases and machine learn the data.</p><p>The crystal structure is optimized in a first optional step (a). This step is optional to ensure that structures of any origin can be used as input for the subsequent COHP analysis. Then, the wave functions are calculated based on DFT (b). In this step, special requisites for the DFT calculation are considered. Due to the required projection from plane waves to atomic orbitals, as many bands must be included in the DFT calculation as are needed for the projection to atomic orbitals. The number of bands actually depends on the structure at hand and on the basis chosen for the projection in LOBSTER. COHPs can only be computed if this number is correctly set by sheer mathematics. In Section 2 of SI, it is demonstrated how the minimum number of bands is computed for GaN. This detail also illustrates the need for automation, since the input in this earlier step depends on an input that is chosen in a later step of the computational procedure, and VASP and LOBSTER runs are highly interdependent. This is also one of the most labor-intensive steps in the entire bonding analysis workflow and requires a lot of training for new users. Our workflow ensures that users always select the correct minimum number of bands according to the proposed basis. Currently, two basis sets (Koga and pbeVASPfit2015) can be used for the entire periodic table but only pbeVASPfit2015 also includes additional orbitals relevant for the solid state (e.g., 2p in metallic Be [17] ). Some more trivial settings for the DFT run must also be chosen. For example, DFT runs must be performed without most of the symmetries to be compatible with LOBSTER (VASP input: ISYM=0) as only time-reversal symmetry is implemented in LOBSTER so far. The next step is to test several atomic orbital bases for the projection of the wave function in LOBSTER (c). Some LOBSTER beginners trained in plane waves but lacking atomic-orbital knowledge may overlook the need to thoroughly test the basis for the projection. Here this is done systematically and all basis functions possible in LOBSTER are exploited. We start with a minimal basis consisting only of occupied valence orbitals in the atomic ground state of each atom. This is determined using the pseudopotential file of VASP. We then systematically add unoccupied orbitals for each element and test all possible combinations. This test is further constrained by the basis functions currently available in LOBSTER. Particular attention has been paid to the handling of the large wavefunction files that are generated as they can become a burden on storage on typical supercomputers, especially for highthroughput studies. They can be automatically deleted. The most important outputs of the VASP and LOBSTER runs will be stored in a MongoDB database (d). For LOBSTER, this is the information from the lobsterout file. Additionally, other output files can be stored in the database as well. In the last step (e), an automatic analysis of the LOBSTER calculations typically with the lowest absolute charge spilling is performed. This charge spilling indicates how much charge was lost in the projection from the plane-wave wave function to the atomic-orbital wave function and it is one of the most important quality criteria for the projection; in principle, the transformation is exact but small deviations may occur, just like in DFT. If used in the "cation-anion" mode, the automated analysis first detects cations and anions based on Mulliken charges calculated using LOBSTER [40] and then identifies the strongest cation-anion bond within the crystal structure based on the integrated Crystal Orbital Hamiltonian populations (ICOHPs). For the test systems shown below, only bonds with a strength of at least 10% of the strongest cation-anion bond were considered. This cutoff can be adjusted. In case all bonds (not only cation-anion bonds) are considered in the analysis, the ICOHP cutoff will typically be set to 10% of the strongest bond. To implement this workflow and simplify the bonding analysis, we extended and created several Python packages. Two new workflows are available in the atomate [28] workflow library (see Figure S1 in the supporting information). The first one is the workflow above. The second one relies only on a predefined basis. This workflow should only be used if the quality of the predefined basis has already been tested. The concrete DFT workflows are implemented in atomate using the workflow manager fireworks, the error handler custodian, and the pymatgen library for material analysis. With atomate and the other underlying Python packages, it is possible to perform complex calculations mostly with VASP (e.g., band structure calculations) on high-performance computers. All input files are generated automatically. To handle the input and output files, we have implemented classes for almost every input and output file of LOBSTER in pymatgen. Plotting classes for COHPs, COOPs, COBIs are now available and there is an interface to the existing plotting classes for densities of states and fatbands in pymatgen. New classes have been developed to allow for the automatic analysis of local environments based on orbital-based bonding analysis. They connect to available tools for determining and analyzing atomic environments in pymatgen (ChemEnv, [36] localenv [37] ). The bonds between atoms are determined based on the magnitude of the ICOHP values. Comparisons with ideal polyhedra are then made based on, for example, ChemEnv and its implementation of continuous symmetry measures describing a distance to a shape [36] . In combination with localenv [37] , structure graph objects and further information about the neighboring atoms can be retrieved. In order to perform hundreds to thousands of calculations in parallel, error handling and validation are critical for highthroughput studies. We have therefore also added new implementations for validators in custodian, which must be successfully passed. For example, this system checks whether the main LOBSTER output file (lobsterout) has been created and whether enough bands have been selected in the DFT run. Also, a validator was added to check if the charge spilling is reasonably low (<5%). We use the latter validator only when exactly one set of basis functions is used for the LOBSTER calculations. If different basis functions are tested in the workflow, this check is not performed. Also, a job class has been introduced to run LOBSTER jobs with custodian. Two tutorials showing how to use these classes and functions are available online . [41,42] Parts of this implementation have been used in other publications, e.g., to understand defects and their stability in photovoltaic materials or to benchmark a new implementation in the LOBSTER program. [16,43] In the latter study, [16] we identified the computational settings used in this study that lead to very well converged LOBSTER results. These and other details on the DFT and LOBSTER settings can be found in the methods section of this paper. The tools for automatically analyzing LOBSTER outputs are implemented in the lobsterpy package. An overview of the package can be found in Figure S1 in the supporting information. The Analysis class takes care of the automatic analysis of LOBSTER output and is the starting point for another class called Description, which will provide automated text and tools for automatic plotting. Detailed examples of COHP data analysis are attached to the code repository corresponding to this paper. By default, only cation-anion bonds are included in the analysis. Especially for large compounds such as Yb14Mn1Sb11, the automatic output otherwise provides too much information for the user, since weak cation-cation interactions are also present. An analysis of all bonds can be enforced since there are compounds where cation-cation interactions are responsible for the stability of a material and the particular structure (e.g., Ge4Se3Te [44] ). For very large coordination environments (coordination number >13), only the coordination number is determined, but not the environment, which is based on a limitation of ChemEnv. ChemEnv only has reference polyhedra for coordination numbers ≤13.</p><p>Lobsterpy also offers a command-line interface that connects and extends the plotting tools available in pymatgen. The user can provide customized style sheets for the plots. This is also possible for automatic COHP plots. The automatic description can also be printed to the screen.</p><!><p>The Lobster workflow in atomate was run for several test systems (NaCl, GaN, MTaO2N (M = Ba, Ca, Sr), Yb14Sb1Mn11). First, we will discuss the resulting projections. For most systems, only one basis is available (see Table 1) since their constituting ions are in noble gas configuration (full octet).  1 shows the Madelung energies per formula unit of GaN. For the COHP plots, a Gaussian smearing was used as the VASP version determined the band gap at the bottom of the band gap with the tetrahedron smearing but not the Gaussian smearing. We expect that this will be fixed again in future versions of VASP. Furthermore, we manually adapted the energy range here and restarted the Lobster run from the projection step. Very often, only states close to the Fermi energy EFermi are relevant for our bonding analysis and we therefore have a different standard setting for the energy range in our automatic workflow.</p><p>For Yb14Sb1Mn11, we will test the influence of the 5d orbital for Yb. This orbital would be unoccupied in a Yb atom according to the corresponding pseudopotential file. However, its use reduces the absolute charge spilling by more than 1%. Therefore, we will use the larger basis in the following. For all systems tested here, the charge spilling is less than 2%, indicating a very good projection. Note that a subsequent Löwdin orthogonalization of the local basis automatically carried out in LOBSTER assures that the entire Hilbert space is recovered, so no electron density gets lost. We assume that this automatic checking of the basis will indeed ensure that users test the projections, rather than simply choosing a basis that gives reasonable spilling under a certain limit without further checking. We also hope that high-throughput studies will allow us to develop further criteria for the quality of a basis besides the absolute charge spilling. We now discuss the results of the automatic analysis of the output files.</p><!><p>We start with the analysis of two rather simple binary systems: GaN and NaCl. We usually describe chemical bonds in crystalline materials in terms of their metallic, covalent, and ionic contributions. The Van Arkel-Ketelaar triangle represents the different contributions for materials and allows us to map materials according to the nature of their bonds. [1,[45][46][47] COHP provides a way to quantify covalent contributions. Madelung energies are commonly used to quantify ionic contributions. LOBSTER calculates these based on quantum-chemical Mulliken or Löwdin charges rather than just based on formal oxidation states. In the following, we will show that such information about the bonding situation in crystals can be easily extracted using our automation. We will also use this information on covalency and ionicity to discuss the stability of the different GaN and NaCl phases.</p><p>We start with the four most stable GaN phases available in the Materials Project database [29] (see Table 2 and Figure 3). These four phases have been investigated by experimental and theoretical studies. [48] GaN is a well-known semiconducting compound for optoelectronic applications and is therefore expected to show significant covalent interactions. [49] Based on COHPs, our code correctly recognizes all coordination environments. This is consistent with geometric evaluations based on ChemEnv and its standard settings. ChemEnv compares coordination environments to ideal shapes. [36] We now proceed to discuss the stability of the compounds and what is effectively causing this. GaN in the wurtzite and zincblende structures is predicted to be more stable than GaN in NaCl structure type based on the DFT total energies, which is consistent with our expectations. GaN in the BN structure type lies between these extremes. We expected that the stabilization of the zincblende and wurtzite structure type compared to the NaCl structure type is due to a higher covalency as GaN is semiconducting. [1] Indeed, we find that the more stable compounds crystallizing in the wurtzite and zincblende types have more covalent interactions than GaN crystallizing in the NaCl structure type. This is indicated by the sum of Ga-N ICOHPs per Ga atom (Table 2 and Figure 3) and it is also connected to the Ga-N distances in the compounds. Ga in wurtzite and zincblende structure type show much stronger COHPs, especially for lowerlying states (mostly with contributions of N (2s)). This is, of course, directly correlated to the smaller distances in these compounds as, for example, the more localized 2s orbitals (around -15 eV below the Fermi energy EFermi) cannot overlap that easily with the Ga orbitals anymore. The dispersion of the bands seems to be very similar in all four compounds, which is a bit surprising. The total energy of GaN in the BN structure, however, cannot be explained that easily based on covalent contributions alone. The calculated Madelung energies suggest that there may be electrostatic stabilization that could explain this energetic ordering. GaN in the BN structure exhibits the most negative Madelung energy per formula unit based on Mulliken charges of all four compounds. The same tendency can be found based on Löwdin charges. As shown in this simple example, correlations of DFT total energies with covalent bond strengths or Madelung energies can now be easily tested on a much larger scale to potentially understand phase transitions based on bonding properties and electrostatics.</p><p>We now turn to the analysis of simple but very ionic systems -two NaCl phases. Here we do not expect covalency to be the driving force for the stability of these compounds. Nevertheless, it is possible to identify the correct coordination environments in these phases based on the bond energy given by COHP (Fig. 4).</p><p>For the binary ground-state structure of NaCl (mp-22862, Figure 4 a), the code correctly recognizes the octahedral coordination environment for Na + and summarizes the mean values for the integrated COHPs. For the high-pressure phase of NaCl (mp-22851, Figure 4 b), the cubic environment is correctly recognized.</p><p>The differences from the ground-state structure of NaCl are visible in the sum of ICOHP values (ICOHP_sum) and the mean value of ICOHPs per bond (ICOHP_mean). The ground-state structure of NaCl has a more negative value (i.e., it has stronger covalent bonds). However, this value should be taken with caution because NaCl is a comparatively ionic compound and exhibits very small ICOHP values in general. In this case, the electrostatics and calculated Madelung energies should lead to a much better understanding of the compound stability: Indeed, the ground-state structure of NaCl also shows a more negative Madelung energy than the high-pressure phase by more than 0.3 eV per cell. Additionally, the automatic output dictionaries for both phases are shown. Note that the automatic output provides a counting scheme for the atoms that starts at 1 (similar to LOBSTER outputs), but that the digits for sites within the dictionaries start at 0 (standard in Python).</p><!><p>Now, two more advanced examples will be discussed. We first focus on the oxynitride compounds CaTaO2N, SrTaO2N, and BaTaO2N, which crystallize in a perovskite-related structure.</p><p>The anion order in these compounds have been under some debate in the past, and possible driving forces for certain anion orders have been discussed. [50][51][52] We now automatically analyze the importance of certain bonds for the anion order in these compounds (Figure 5). For this purpose, we vary the anion order of the systems in a supercell with up to 2 formula units. This results in seven different structural models per compound with different anion orders. The energetically least and most favorable models are shown in Figure 5 a and b. Based on our automated bonding analysis, we can then easily correlate the total energy of the systems with the strongest covalent interaction (Ta-N bond) in the systems (Figure 5 c). We also see an anti-correlation with bond energies of the Ta-O bonds because the more covalent Ta-N bonds strengthen at the expense of the Ta-O bonds. The covalency of the Ta-N interactions is therefore crucial for the anionic ordering of these oxynitrides.</p><p>The Ca/Sr/Ba-O bonds are too weak to be detected with our automated analysis and the corresponding relative ICOHP cutoff (10% of the absolute largest ICOHP). The ICOHPs help to focus on the most important covalent interactions in our system, which could also be an advantage for determining coordination environments compared to purely geometric determinations of the latter. Thanks to our automation tools, the whole bonding analysis process can be performed with only a few lines of Python code, as shown in the repository accompanying this paper. This correlation between the strongest covalent interaction and the total energy of the system has already been observed for the oxynitride CaTaO 2N [50,51] . This correlation has now been confirmed for SrTaO2N and BaTaO2N, which was to be expected given the chemical similarity of the compounds. One now arrives at such an analysis almost fully automatically. As a final example, we show the analysis of a very large and complex structure, namely that of the well-known thermoelectric material Yb14Mn1Sb11. [53] It has 104 atoms in the unit cell, and such a bonding analysis by hand would require a lot of manual work. Up to a distance limit of 6 Å (maximum considered distance in our current implementation of the workflow), more than 1400 bonds are found in this structure. The band structure of the compound in the nonmagnetic state has been recently studied. [54] Similarly, an analysis of the band structure of Yb14Mg1Sb11 crystallizing in the same type of structure, has been carried out. [55] The magnetism in Yb14MnSb11 has also already been studied both experimentally and theoretically. [56] From previous studies, we expect an influence of the Mn-Sb bonds on the magnetism.</p><p>To investigate the influence of chemical bonding on magnetism in this compound, we decided to compare COHPs from nonmagnetic and ferromagnetic runs. This type of investigation was inspired by a study by Landrum and Dronskowski on magnetism and driving forces for itinerant magnetism based on chemical bonding. [15] We then performed two static LOBSTER runs that differ in the magnetic model after structural optimizations in the ferromagnetic setting. We introduce magnetic moments on Mn only. Similar to Perez et al. [54] in the nonmagnetic state, Yb, Mn, and Sb states play an important role in the valence band of Yb14Mn1Sb11. We find that both Yb-Sb and Mn-Sb bonds are relevant cationanion bonds in the structures, the Yb-Sb bonds being on average stronger than the Mn-Sb bonds. The bonds in the "Sb3 7-" are not considered in the standard automatic analysis of LOBSTER output files. By default, this analysis focuses only on cation-anion bonds to obtain a readable result even for very large structures such as Yb14Mn1Sb11. However, it is possible to enforce their analysis. The automatic text output and the structure for the ferromagnetic setting are shown in Fig. 6.  [57] and automatic analysis of the conventional cell for CaTaO2N. The code is able to automatically determine relevant bonds, coordination environments and average bond strengths (average ICOHPs); note that more negative ICOHPs indicate stronger bonds. In this case, only the Ta-O and Ta-N bonds were detected. Potential Ca/Sr/Ba-O and Ca/Sr/Ba-N bonds are much less covalent and play a less important direct role for the material. This information can then be easily correlated with other material properties (e.g., total energies) to understand these properties based on bond strengths. b) Structure of the most stable Ca/Sr/Br TaNO2 model including up to 2 formula units.. c) A correlation plot of the most important ICOHPs for this compound and the total energy of the systems with different anion orders is shown. The correlation coefficient r for the strongest covalent interaction (Ta-N) is always higher than 0.95, indicating a significant correlation. In the ferromagnetic setting, the octahedral environments are discovered for three of four crystallographically nonequivalent Yb atoms. The other is determined to be a pentagonal bipyramid. In contrast, the purely geometric tool ChemEnv, [36] based on Voronoian analysis with distance and solid angle cutoffs and determination of continuous symmetry measures, determines octahedral environments for all four. SimplestChemEnvStrategy and a standard relative solid angle cutoff of 0.3 were used. This lower cutoff is relative to the largest solid angle in the Voronoi construction. One of the identified octahedral environments is, however, highly distorted, as indicated by a continuous symmetry measure greater than 7. This is the environment that our orbitalbased analysis determines to be a pentagonal bipyramid. If we change the angle cutoff to 0.2, the geometric algorithm also arrives at the pentagonal bipyramidal environment. The determination of the tetrahedral coordination environment for Mn agrees well with the geometric-based determination. Overall, the agreement between the geometric and orbital-based determination of the coordination environments is very good.</p><p>Upon analyzing the COHP diagrams with and without spin polarization for Yb 14Mn1Sb11 (Figure 6b), we indeed find an influence of the Mn-Sb bonds on the magnetism as it has been seen before [56] : strong antibonding Mn-Sb interactions are discovered below the Fermi level for the calculation without spin polarization. When spin polarization is turned on, the antibonding interactions in one spin channel disappear completely and this spin channel now has fully bonding Mn-Sb bonds. In the other spin channel, small antibonding Mn-Sb interactions remain around the Fermi level. This finding agrees well with the antibonding fingerprint for other ferromagnetic compounds such as the transition metals, FeNi3, FePd3, and MbSb, demonstrated by Landrum and Dronskowski. [15] In contrast to the results of Landrum and Dronskowski, we do not detect stronger Mn-Sb bonds for the ferromagnetic setting, however. These bonds are slightly weaker than in the non-magnetic setting. In our gedankenexperiment, however, we find an "oxidation" of Mn when spin polarization in considered. The Mulliken charges of Mn change dramatically from 0.69 (nonmagnetic) to 0.91 (ferromagnetic). This makes the Madelung energy much more negative with a change from -104.21 to -118.82 eV. This could indicate additional electrostatic stabilization for this compound in the ferromagnetic setting, rather than the typical covalent stabilization observed for other compounds such as α-Fe. These differences to the situation in α-Fe with homopolar Fe-Fe bonds are not surprising as we now observe this change in a heteropolar bond where Sb is more electronegative than Mn and thereby captures some of the electrons of Mn when spin-polarization is switched on.</p><!><p>This work simplifies bonding analysis by providing workflow tools and tools for automated analysis. Prior to this study, chemical bonding analysis based on orbitals was a task with many manual steps and therefore error prone. We have demonstrated our automation here using several test systems, including the study of the anionic order of several oxynitrides and the complex crystal structure of Yb 14Mn1Sb11. Thanks to automation, these calculations can now be performed in a high-throughput manner.</p><p>In addition, we have shown how the tools allow correlating important bonding properties with other material properties. We expect that these tools will facilitate the search for new descriptors for machine learning of material properties. Bader charges [58] or other information based on electronic structure theory have already helped in such studies. In addition, we hope that the tools will help to provide further chemical understanding of materials.</p><p>The study also illustrates how much work is needed to automate just one task in computational chemistry. Future work will aim to include orbital-by-orbital analysis of bonds, various indicators of covalency (ICOOP, ICOBI), and k-dependent covalency in our automated analysis, which could be particularly useful for compounds with very steep bands. In addition, pre-convergence steps with a smaller number of bands could be included in our LOBSTER workflow to speed up calculations and convergence of calculations, which is currently a problem for compounds requiring a very large number of bands (such as Yb14Mn1Sb11).</p><p>Interfaces to other DFT codes could be developed.</p><p>For all examples shown, we provide the analysis scripts in our Github repository (https://github.com/JaGeo/LobsterAutomation, link to zenodo.org will be provided in the final publication).</p><!><p>For all materials considered in this study, DFT calculations based on the PAW method [59,60] and the PBE approximation of the exchange correlation functional [61] were performed. For optimization, we relied on many of the default settings of the atomate package imported via pymatgen from the MPRelaxSet. However, we chose more precise optimization settings. We use an energy difference criterion for the convergence of the electronic structure of 10 -6 eV, and the structure is considered relaxed if the energy difference between two consecutive steps is below a setting of 10 -5 eV. We relied on the standard cutoff for the plane-wave energy of 520 eV. In contrast to the MPRelaxSet in pymatgen, we used the VASP pseudopotentials with version number 5.4, with the additional difference that s-electrons are included in the valence of W. We used 6000 k-points per reciprocal atom to ensure that the energies for the defect phases of the oxynitrides converged well. We also chose a Gaussian smearing for the optimization. When calculating the wave function, the energy is converged until it is less than 10 -6 eV. This stricter criterion, as opposed to the default setting in the materials project, was necessary to ensure that the wave functions were of sufficient quality. Otherwise, we would sometimes have obtained different charge spillings for the two spin channels of actually not magnetic compounds, which is contrary to expectation and indicates problems with convergence of the wave function. In addition, we considered time-reversal symmetry in the VASP calculation and used 6000 k points per reciprocal atom to ensure convergence of the LOBSTER results.</p><p>The number of bands is adjusted according to the largest basis available in LOBSTER for the compound in question. The other parameters are the established defaults for the Materials Project PBE calculations as implemented in MPRelaxSet. The entire implementation of a generator for the input can be found in pymatgen and in pymatgen.io.vasp.sets.LobsterSet. These settings were previously determined to produce well-matched results when testing a new implementation in LOBSTER for several hundred LOBSTER calculations. [16] The LOBSTER run itself was then based on the pbeVASPfit2015 basis and fits the atomic basis to the selected POTCAR files. All standard outputs of a LOBSTER run (i.e., information on Crystal Orbital Hamilton populations, Crystal Orbital Overlap Populations, Crystal Orbital Bond Indices, projected density of states, Mulliken charges, Löwdin charges, etc.) were calculated.</p><p>To generate the different anion configurations of the perovskite structures, we again relied on pymatgen, enumlib [62] and also on symmetry identification based on spglib [63] .</p><p>The following program versions have been used: pymatgen 2022.2.1, atomate 1.0.3, LobsterPy 0.1.0, Lobster 4.1.0, VASP 6.2.1.</p>