Crystal orbital overlap population based on all-electron ab initio simulation with numeric atom-centered orbitals and its application to chemical-bonding analysis in Li-intercalated layered materials

Crystal orbital overlap population (COOP) is one of the effective tools for chemical-bonding analysis, and thus it has been utilized in the materials development and characterization. In this study, we developed a code to perform the COOP-based chemical-bonding analysis based on the wave function obtained from a first principles all-electron calculation with numeric atom-centered orbitals. The chemical-bonding analysis using the developed code was demonstrated for F2, Si, CaC6, and metals including Ti and Nb. Furthermore, we applied the method to analyze the chemical-bonding changes associated with a Li intercalation in three representative layered materials: graphite, MoS2, and ZrNCl, because of their great industrial importance, particularly for the applications in battery and superconducting materials. The COOP analysis provided some insights for understanding the intercalation mechanism and the stability of the intercalated materials from a chemical-bonding viewpoint.


Introduction
Chemical-bonding among atoms composing materials is closely linked to their stability, reactivity, and electronic properties, and it has been investigated using simulations in a first principles manner.
In other words, we can delve more deeply into the mechanisms underlying the emergence of materials properties by extracting information regarding chemical bonds.For the chemical-bonding analysis, some theoretical methods have been utilized, such as density of states (DOS), partial/projected DOS, and direct visualization of the wave function.One of the powerful methods in computational chemicalbonding analysis is Crystal Orbital Overlap Population (COOP) analysis [1], where orbital overlap is quantified to understand bonding and anti-bonding characteristics of chemical-bonding.Furthermore, similar concepts such as Crystal Orbital Hamiltonian Population (COHP) [2] have been devised and utilized in this context.Over the years, various DFT calculation codes have been created, whereas the development of codes for post-process COOP/COHP calculations has still been limited.The most straightforward approach to perform COOP/COHP computations is to utilizing wave function obtained from DFT calculations using localized basis sets.Computer programs for such COOP/COHP calculations have been implemented in DFT codes based on Linear Combination of Atomic Orbitals (LCAO) method like SIESTA [3,4] where atomic orbitals are employed as basis sets.Furthermore, a program named Local Orbital Basis Suite Towards Electronic-Structure Reconstruction (LOBSTER) [5][6][7] enables chemical-bonding analysis based on the results from plane-wave basis DFT simulations using The Vienna Ab initio Simulation Package (VASP) [8,9], Quantum Espresso [10], and ABINIT [11], by projecting plane-wave basis wave-functions onto localized basis sets.However, these methods have been developed based on non-all electron method (such as pseudopotential method) and/or delocalized basis method (such as plane-wave basis method).If we can calculate COOP based on the all-electron DFT calculation with atom-centered orbitals, a "direct" chemical-bonding analysis without any projection step including core and semi-core sates become available.
In this research, to further facilitate the chemical-bonding analysis of materials and molecules, we have newly implemented a code to calculate the COOP for arbitrary orbital pairs and arbitrary atomic pairs based on the wave-function obtained by FHI-aims [12], which is an all-electron DFT code employing numeric atom-centered orbitals as basis functions.Using the developed program, we have demonstrated the chemical-bonding analysis on simple molecular and solid-state system, F2 and diamond type Si, and confirmed the validity and effectiveness.
Further, as an application of COOP analysis in materials research, we have performed a chemicalbonding analysis in Li intercalated layered materials.The phenomenon of Li intercalation into layered materials is not only utilized in battery and superconducting materials but also in modulating electrical and optical properties of materials [13][14][15].It is an industrially significant phenomenon and gaining a deeper understanding of the mechanism behind Li intercalation is of great importance.In this study, we have analyzed the changes in chemical-bonding state associated with Li intercalation for representative layered materials ranging from one-elemental to ternary systems: graphite, MoS2, and ZrNCl, to gain new insights into the phenomenon.We selected these systems because understanding the changes in electronic structure associated with the Li insertion is essential for developing industrially important functions in materials such as batteries and superconductors.Through the COOP based chemical-bonding analysis, we show the effectiveness of our codes.

Definition of COOP
The COOP, which stands for "Crystal Orbital Overlap Population", is based on the extension of the Mulliken population analysis to assign electrons to atoms, to assign electrons to chemical bonds in molecular orbital method [16,17].COOP analysis is one of the useful and representative methods for qualitatively and quantitatively analyzing chemical-bonding in materials.
The Kohn-Sham equation, which needs to be solved self-consistently in DFT calculations, for one electron wave function | , � in periodic system is: where the electronic Hamiltonian ℎ � KS is written as: Here,  ̂ ,  �  ,  �  ,  �  are the Kohn-Sham effective single-particle kinetic energy, the external potential, the electrostatic potential of the electron density, and the exchange-correlation potential, respectively.In addition,  ,  , and   () are wave vector, band index, and the corresponding eigenvalue.In non-periodic system, the Kohn-Sham equations and their solutions are not dependent of k.When | , � is the basis function, the solution to the Khon-Sham equation can be written as: where  is the index of basis function and  , () is the corresponding expansion coefficient.In Mulliken's approach, the overlap population  , () *  , () , , which is calculated using the overlap matrix element  , = � , � , �, is defined to determine the electrons associated with the  -th and -th orbital.The DOS weighted by the overlap population associated with -th orbital is called partial density of states (PDOS): , which enables us to obtain the contribution of -th orbital to the DOS.In practical calculation, real part of the overlap population is used for the weighting.By multiplying the occupancy number and performing integration, it is possible to determine the number of electrons belonging to -th orbital, which is defined as the Mulliken charge belonging to  -th orbital.The sum of Mulliken charge belonging to an atom corresponds to the effective charge of the atom, and the subtraction of atomic number provides the net charge of the atom.
Similarly, by considering only the pair of orbitals  and , and weighting the DOS with the overlap population, we can obtain: The sign of COOP is related to the type of interaction between orbital  and , where a positive COOP corresponds to a bonding interaction and a negative COOP corresponds to an antibonding interaction.
In other words, COOP diagram provides a visual representation of the bonding and anti-bonding characteristics of orbital interactions, allowing us to understand the nature of chemical-bonding in materials.The value obtained by integrating the COOP multiplied by the occupancy is called the bond overlap population (BOP), which corresponds to the number of electrons shared between orbital  and , and can be considered a quantity that is correlated with the strength of covalent bonds.

Computational details
Prior to the COOP analysis, the local structure of the materials was commonly optimized by a firstprinciples calculation using the plane wave basis Projector Augmented Wave (PAW) method [18], as implemented in VASP package [8,9] because some of their structures are not experimentally reported.
The structural optimizations were carried out using the Generalized Gradient Approximation by Perdew, Burke, and Ernzerhof (GGA-PBE) [19] for the treatment of exchange-correlation interactions in F2 molecule and diamond Si, while rev-vdW-DF2 [20][21][22][23] by Hamada was adopted to account for the van der Waals interactions in the case of graphite, MoS2, ZrNCl, and their structures with Li insertion.The cutoff energy for the plane-wave expansion was set to 650 eV.For k-point sampling, a k-mesh with a spacing of 0.25 Å -1 was used, with exception of F2, for which only a single gamma point was employed.Ionic relaxation was carried out until the forces acting on each atom were reduced to less than 0.01 eV/Å.
To investigate the correlation between the chemical-bonding and intercalation behaviors, intercalation energy was calculated using the following equation: Here,  [Li] represents the total energy of an isolated Li atom, and [LM] represents that of a host layered material (LM), graphite, MoS2, and ZrNCl in this case.[LM with Li] is the total energy of the Li-intercalated layered materials.
After performing the structural optimizations, a chemical bonding analysis was conducted using FHIaims [12].The electron-electron interactions were treated using GGA-PBE, and the calculation were performed at twice the k-point density of the structural optimization to obtain converged DOS/PDOS/COOP diagrams.After the self-consistent field simulation, the COOP (which is OP for molecule) was calculated based on the overlap population matrix elements  , (  for molecule) and expansion coefficients  , () ( , for molecule).A minimal basis set was employed to enhance interpretability, whereas a light basis set was utilized for CaC6, graphite, MoS2, and ZrNCl systems, allowing for analysis up to unoccupied orbitals while still maintaining interpretability.

F2 molecule
First, to demonstrate the use of COOP in chemical-bonding analysis, we consider the example of a simple molecule, F2.Since molecular system is not "crystal", we call it overlap population (OP) diagram here.In addition, the interaction between 2  orbitals result in the formation of a bonding  3 and antibonding  3 * orbitals, while the interaction between 2  (2  ) results in the formation of bonding  and anti-bonding  * orbitals.Based on the definition of OP, we expect that the OP values to be positive for bonding orbitals such as  and , and negative for anti-bonding orbitals such as  * and  * .

Diamond type Si
The second example is diamond-type Si, which is one of the representative covalent solid-state materials.Figure 2 shows the calculated PDOS and COOP diagrams.The results obtained using FHIaims are presented in the leftmost column (a).To test the validity of the results, additional calculations based on the LCAO method were performed using SIESTA [3,4] and the plane-wave basis method using VASP [8,9], as shown in columns (b) and (c), respectively.SIESTA is a first-principles calculation code that utilizes numerical atomic orbitals, employing the pseudopotential method.The COOP of the identical Si crystal was obtained using a minimal basis set, with the k-grid cutoff set at 30Å.The VASP simulation was performed with a cutoff energy of 650 eV and a k-mesh of 15×15× 15, and the PDOS/COOP was obtained through post-processing calculations using the LOBSTER [5][6][7].
The first row of Fig. 2 represents the l-resolved PDOS for a Si atom in Si crystal, and the figures on the second row shows the calculated COOP diagram regarding Si-Si interactions.In the third row, we show the COOP diagrams regarding 3s-3s, 3s-3p, and 3p-3p interactions.By comparing these results, it can be confirmed that the overall shapes are in good agreement across all methods.In particular, FHI-aims and SIESTA, which utilizes numeric atom-centered basis function, show very similar PDOS and COOP in both valence and conduction bands.From the COOP diagrams, we can confirm that the valence band is formed by Si-Si bonding orbitals, while the conduction band is formed by strong Si-Si anti-bonding orbitals.
Furthermore, l-resolved COOP can elucidate the detailed chemical-bonding characteristics of Si.In the valence band, there are three peaks, A, B, and C. Peak A, which is mainly composed of Si-3s components, is primarily attributed to the 3s-3s bonding interactions and small components of the 3s-3p bonding interactions.Peak B, which is composed of both 3s and 3p components, is attributed to the 3s-3p bonding interaction and 3s-3s antibonding interaction, whereas the valence band maximum peak C to the bonding interaction between 3p-3p orbitals.
In this way, also in solid state material, the COOP diagram can be used to visually comprehend the characteristics of the orbitals in the bands, and thereby deepen the understanding of chemical-bonding states.Although the different simulation codes show consistent PDOS/COOP in the diamond-type silicon case, the same PDOS/COOP diagrams are not necessarily guaranteed in some cases.For instance, Fig. 3 shows the results for CaC6.This figure includes the total DOS, the sum of partial/projected DOS, species-specific PDOS, and COOP diagrams.Figure 3 (a) details the results computed using FHI-aims, Fig. 3 (b) outlines results from VASP+LOBSTER, and Fig. 3 (c) shows the results derived solely from VASP.Generally, when calculations are performed using a plane-wave basis set and wave functions are projected onto atoms, the sum of projected DOS does not need to identical to the total DOS, as illustrated in Fig. 3 (c).LOBSTER uses a refined projection scheme to mitigate this issue; however, when a suitable local basis is missing within the provided basis set, wave functions may be mistakenly projected onto incorrect atoms (Fig. 3 (b)).
Specifically, for CaC6, the first sharp peak A in the conduction band is originating from the 3d component of Ca atom.If these orbitals are not property included, the projection could inaccurately occur onto C atoms instead as shown in the second row of Fig. 3(b), leading to the inaccurate PDOS and COOP diagrams where anti-bonding interactions between C and Ca is not properly evaluated (the third row of Fig. 3 (b)).In our method, which employs numeric atom-centered basis function for electronic structure relaxation, it is possible to include any orbital in the electronic structure calculation, allowing the results to be directly utilized in the PDOS and COOP analysis.This capability enables the flexible chemical-bonding analysis across a wide range of materials.

Applications: Chemical-bonding analysis in Li-intercalated Layered Materials
As an example of the application of COOP analysis, in this study, we analyzed the changes in chemical bonding states associated with the intercalation of Li atoms in layered materials.We selected three representative layered materials, namely graphite, MoS2, and ZrNCl, as shown in Fig. 3, ranging from one-elemental system to ternary systems.

Graphite
Graphite consists of stacked graphene layers held together by van der Waals interactions.Due to the relatively weak bonding between graphene layers, it easily allows the insertion of atoms or molecules into the layers [25].The most typical application of atomic intercalation in graphite is its use as an anode material in Li-ion batteries.During the charging process, Li intercalates into the interlayer spaces of graphite, forming a graphite intercalation compound known as LiC6.LiC6 can stably exist, and its origin has been computationally investigated so far [26][27][28].However, the specific changes in chemical bonding states associated with the intercalation process have been rarely investigated.
In this study, calculations were performed for three models with different Li concentrations: C, Li0.021C (LiC48), and Li0.167C (LiC6).We have focused only on the AA stacking structure, where both the layers that contain Li and those without Li exhibit AA-stacked structure, as it is known that the stacking structure of graphite changes from AB stacking to AA stacking upon Li intercalation.In Li0.021C, Li atoms are intercalated every two graphene layers, which is referred to as stage-2 structure.
In the case of Li0.167C, Li atoms are intercalated between all layers as shown in Fig. 4 (a).Li and C, with charge density maps serving as the evidence for this argument [26].In this study, the COOP analysis successfully quantitated the bonding interactions between Li and C, and revealed that the covalent bonds gradually strengthen as the concentration of Li increases.On the other hand, the COOP diagram also reveals that the electrons start to occupy the C-C anti-bonding orbitals (second row of Fig. 6), which likely weaken the covalent bonds between them.In fact, the BOP values of C-C interaction decrease from 0.500e in graphite to 0.460e in Li0.167C.The intercalation energy increases from -2.213 eV in Li0.021C to -1.887 eV in Li0.167C, and the increase in Li content contributes to destabilization, which is possibly correlated with the destabilization of C-C bonds in terms covalent bonding.
To summarize, from the above results, it is considered that when Li is intercalated in graphite, the charge is transferred from Li to C from the viewpoint of Mulliken charge, and an ionic interaction occurs.In addition, COOP analysis revealed that the intercalation of Li induces covalent bonding interactions between Li and C, while anti-bonding interactions occur between C atoms, resulting in the destabilization of the structure.to be per C atom.

2H-MoS2
MoS2 is a representative two-dimensional compound of binary systems and belongs to the class of transition metal dichalcogenides (TMDs).Similar to graphite, Li can be inserted between the layers of MoS2, and the phenomenon of Li intercalation is utilized for the applications in battery materials, modulation of electrical properties, and exfoliation of MoS2 layers [13].It is known that MoS2 undergoes successive structural phase transition from 2H to 1T, and then to 1T' phase [29].In this research, we placed emphasis on the initial phase, 2H phase, to investigate the interaction between Li and the MoS2 layer in detail.Specifically, we focused on the region with relatively low Li concentrations: MoS2, Li0.056MoS2, and Li0.250MoS2, and examined the changes in the chemicalbonding states associated with the Li interactions.In Li0.056MoS2, Li is inserted every two layers, while in Li0.250MoS2, which has the highest Li concentration, Li is inserted every interlayer, as shown in Fig.In the case of MoS2, it was observed that Li insertion led to the transfer of electrons from Li to S.

(b)
However, even with an increase in Li concentration, the COOP diagrams for both Li-S and Mo-S interactions were not significantly affected.Furthermore, the change in the intercalation energies was also small.This indicates that the changes in the chemical bonding state are closely related to the stability of the intercalated structures.Fermi energy (EF).The total DOS has been scaled to be per MoS2.

ZrNCl
Finally, we have applied the COOP analysis for a ternary layered material, ZrNCl, which can also allow Li intercalation between the ZrNCl layers, leading to the emergence of superconductivity [30,31].
Because it is known that stacking structure changes from SmSI type to YOF type upon Li intercalation [32], we considered YOF type structure with three different Li concentrations: ZrNCl, Li0.042ZrNCl, and Li0.500ZrNCl.In Li0.042ZrNCl, Li is inserted every three layers, while in Li0.500ZrNCl, Li exists in every interlayer as shown in Fig. 4 (c).
Calculated total DOS and atom PDOS are shown in Fig. 9.The figures on the first row represents the total DOS, and the second and third rows represent the PDOS of Zr and Cl bonded to Li, respectively.
In the case of Li0.042ZrNCl, the dilute in-plane configuration of Li results in several symmetrically distinct atom sites for Zr, Cl, and N, and we selectively investigated PDOS of atoms close to intercalated Li.From PDOS in Fig. 8, it can be observed that the electronic states of all elements undergo changes with an increasing Li content, and changes in the peaks occur at the same energy levels in the PDOS of different elements as indicate by the gray solid lines.
From the perspective of Mulliken charges, that of Li does not change significantly even with an increase in Li.The Mulliken charge of Li is +0.417e and +0.436e in Li0.042ZrNCl and Li0.500ZrNCl, respectively.The Mulliken charge of Cl, which is in contact with Li, changes from -0.312e before insertion to -0.434e in Li0.500ZrNCl, and thus Cl tends to acquire a more negative charge.In addition, not only Cl but also Zr undergo changes in Mulliken charge.In ZrNCl, Zr has +1.032e, while it decreases to +0.892e in Li0.500ZrNCl, acquiring a negative charge approximately equivalent to that received by Cl.However, the COOP values between Li and Cl around -9~-7eV decrease with an increasing Li content.
As a result, the BOP values slightly decrease from +0.082e in Li0.042ZrNCl to +0.077e in Li0.500ZrNCl.This is likely be attributed to the high electronegativity of Cl, suggesting that the Li-Cl bonding becomes more ionic by increasing the Li concentration.
As for the COOP between Zr and Cl, it is confirmed that bonding interactions between Zr and Cl decreases, and anti-bonding interactions occur with Li insertion.Consequently, the BOP values between Zr and Cl exhibit a monotonic decrease from +0.130e in ZrNCl, +0.099e in Li0.042ZrNCl, and further to +0.058e in Li0.500ZrNCl.In contrast, there is significant enhancement of the Zr-N bonding interaction at the energy denoted by "p1", where the Zr-Cl anti-bonding interactions are also strengthened.This increase of the Zr-Cl anti-bonding interaction would influence to the intercalation energies.The intercalation energy increases from -2.893 eV in Li0.042ZrNCl to -2.218 eV in Li0.500ZrNCl, which is thought to be correlated with the destabilization of the Zr-Cl bonding.

Conclusion
In this work, we developed a code to perform COOP analysis based on all-electron DFT calculation with numeric atom-centered orbitals using FHI-aims.Then, we conducted the COOP analysis on molecular and solid-state systems: F2, Si, and CaC6, confirming the validity of the present method and demonstrating the chemical-bonding interpretation.Furthermore, we conducted the COOP analysis to investigate the changes in the chemical-bonding associated with a Li intercalation in three representative layered materials: graphite, MoS2, and ZrNCl.
In graphite, it was clarified that Li intercalation leads to bonding interactions between Li and C, while C-C anti-bonding interactions emerge, leading to the destabilization of the host structure.As for 2H-MoS2, while there was charge transfer from Li to S, the change in the chemical-bonding state was small.This small chemical-bonding change was consistent to the small change in the intercalation energy of Li to 2H-MoS2.In ZrNCl, electronic states of all elements were affected by the Li intercalation, weakening the covalent bond between Cl and Zr.On the other hand, it was revealed that new bonding orbitals are formed between Zr and N, effectively capturing the effects of Li intercalation.
In conclusion, the COOP analysis enabled us to accurately understand the change in the characteristics of chemical-bonding associated with Li intercalation and elucidate the origin of stability in Li intercalated materials from a chemical bonding perspective.

Figure 1 (
a) schematically illustrates the bonding and anti-bonding interactions of F2 based on the simple chemical-bonding concept. and  indicate bonding interactions, while the addition of "*" denote antibonding.Based on the simple chemical-bonding concept, the interaction between 2 orbitals of the two F atoms give rise to bonding orbital  2 and anti-bonding orbital  2 * .

Figure 1 (
Figure 1 (b) shows the calculated l-resolved PDOS for a F atom in the F2 molecule, along with the

Figure 1
Figure 1 (a) A schematic diagram depicting the formation of molecular orbitals from atomic orbitals in F2 molecule.Sphere and dumbbell represent the 2s and 2p orbitals, respectively, and their phase was represented by blue and white colors.The white circle shows the electron.(b) Calculated lresolved PDOS and COOP diagram in F2.The vertical dashed line at zero represents the position of Fermi energy, EF.

Figure 2
Figure 2 (a) l-resolved PDOS for single Si atom, COOP diagram for Si-Si interactions, and COOP diagram illustrating the interactions between respective orbitals, calculated using FHI-aims.(b) Results from SIESTA code, and (c) Results from VASP+LOBSTER code.The vertical dashed line corresponds to the position of calculated Fermi level, and vertical solid lines show the positions of peak A, B, and C.

Figure 3
Figure 3 (a) From top to bottom: Total DOS with sum of PDOS, atom-PDOS, and C-Ca COOP diagrams for CaC6 calculated using FHI-aims, (b) VASP+LOBSTER, and (c) VASP.The vertical dashed line at zero represents the position of Fermi energy and the solid line shows the position of peak A.

Figure 5
Figure 5 shows the calculated DOS and PDOS of graphite, Li0.021C, and Li0.167C.The figures on the

Figure 5
Figure 5 Calculated total DOS and atom-PDOS in (a) graphite, (b) Li0.021C, and (c) Li0.167C.The vertical dashed line corresponds to the position of Fermi energy (EF).The total DOS has been scaled to be per C atom.

Figure 6
Figure6shows the results of COOP analysis at each Li concentration.The first, second, and third

Figure 6
Figure 6 Calculated total DOS and COOP curves in (a) graphite, (b) Li0.021C and (c) Li0.167C.The vertical dashed line corresponds to the position of Fermi energy (EF).The total DOS has been scaled

Figure 7
Figure7shows the calculated total DOS and PDOS for each composition.The first, second, and third

Figure 7
Figure 7 Calculated total DOS and atom-PDOS for Mo, S, and Li in (a) MoS2, (b) Li0.056MoS2, and (c) Li0.250MoS2.The vertical dashed line corresponds to the position of Fermi energy (EF).The total DOS has been scaled to be per MoS2.

Figure 8
Figure 8 shows the results of COOP analysis.Total DOS, COOP diagram about Mo-S interactions,

Figure 8
Figure 8 Calculated total DOS and COOP diagrams about Mo-S interactions and Li-S interactions in (a) MoS2, (b) Li0.056MoS2, and (c) Li0.250MoS2.The vertical dashed line corresponds to the position of

Figure 9
Figure 9 Calculated total DOS and atom-PDOS for Zr, Cl, N, and Li for (a) ZrNCl, (b) Li0.042ZrNCl, and (c) Li0.500ZrNCl.The vertical dashed line corresponds to the position of Fermi energy (EF).The total DOS has been scaled to be per ZrNCl.
Based on the aforementioned discussion, the way to accept the Li intercalation by changing the electronic structure of ZrNCl is summarized as follows: Upon the intercalation of Li in ZrNCl, Cl and Zr accept electrons from Li, weakening the covalent bond between Cl and Zr.At the same time, Zr and N are forming a new chemical bonding state.Namely, the elements constituting the layered compounds synchronously change their chemical bonding and accept the Li intercalation, which was first revealed by the COOP analysis.The energy for intercalation increased with the increase in Li content, which is possibly correlated with the destabilization of the bonding between the Zr and Cl from the perspective of covalent bonding.

Figure 10
Figure 10 Calculated total DOS and COOP diagram about Zr-Cl interactions, Zr-N interactions, and Li-Cl interactions for (a) ZrNCl, (b) Li0.042ZrNCl, and (c) Li0.500ZrNCl.Vertical dashed lines correspond to the position of Fermi energy (EF).The total DOS has been scaled to be per ZrNCl.