Noncomparative scaling of aromaticity through electron itinerancy

Aromaticity, a topic of enduring interest in Chemistry, is one such concept which lacks a proper definition and the term is used by most of the chemists with an intuitive notion.¹ Introduced by Kekule in 1865, aromaticity has emerged as an important phenomenon for its intimate relation with chemical reactivity.² The increasing interest on aromaticity developed an urge to attach a quantitative character to it beyond its hitherto available qualitative description. As a consequence, a number of quantitative aromaticity indices have been proposed starting from the famous 4n + 2π-electron rule of Huckel to the most recent Shannon Entropy index.³ At present, the molecular property based criteria for aromaticity can be broadly divided into four basic categories such as structural,⁴ electronic,⁵ magnetic,⁶ and energetic.⁷ Later, electron density based indices such as the para-delocalization index (PDI),⁵ fluctuation index (FLU),⁸ and multi-centre bond index (MCBI)⁹ have been introduced to characterize aromaticity. In a different note, Klein and co-workers correlated aromaticity of molecular benzenoids to a unique type of polynomial,¹⁰ termed as Clar-2-nomial on the basis of Clar’s concept of aromatic sextet.¹¹ Later on, this correlation has been extended for radical benzenoids,¹² and more recently on polyacenes and beyond.¹³ However, most of the indices are being used as indirect measure of aromaticity, as the phenomenon is not directly observable.¹⁴ After a thorough revision of structural, energetic and magnetic indices of aromaticity for many aromatic systems, Cyranski et al. concluded aromaticity as an abstract idea from philosophical point of view.¹⁴ Whereas, from another perspective, the phenomenon of aromaticity can be regarded as multidimensional since it cannot be described unambiguously by any single property based index and all indices do not always give consistent result. Owing to all these issues the quantification of aromaticity still remains a challenging and fundamental problem,¹⁵ and requires to be addressed through a precise phenomenological description.

In spite of different approaches to measure aromaticity, it can be undoubtedly stated that all the aromaticity indicesstem from the delocalization of electrons in a closed loop of favourable symmetry.^1(c),3 In contrast, the unfavourable symmetry of antiaromatic systems leads to localized, rather than delocalized electron structures.¹⁶ In fact, by 1960s, most chemists agreed to define aromatic compounds as cyclic planar systems with delocalised π-electrons.¹⁷ The energy minimization due to this electron delocalization in aromatic systems is often taken as preferred aromaticity index,^15,18 as because the energy criterion of aromaticity governs reactions and much of the chemical behaviour. Hence, in this work, we set our objective to devise a phenomenological model through which the aromaticity can be quantified in terms of the extra stabilization gained due to cyclic delocalization of π-electrons. Identifying this specific part of energy has been an intriguing problem to theoreticians.¹⁹ For the archetypal aromatic systems such as benzene; Pauling and Whealand prescribed to use the energy difference between a hypothetical single valence band (VB) structure and a real multi-valence band structure as an estimate of π-electron delocalization energy (DE).²⁰ They named this energy as Resonance Energy (RE) and applied empirical procedures to obtain the numerical value of RE or in a subtle way – aromatic stabilisation energy (ASE).²¹ Breslow proposed to use the energy difference between the concerned ring and its open chain analogue as a measure of the ASE.²² The overestimated ASE in this process has been rectified by Malrieu and co-workers by considering a double cut in the ring, instead of a single cut.²³ Lenthe and co-workers adopted ab initio VB approach and VB self consistent field method to optimize 1,3,5-cyclohexatriene to compare its energy with benzene for calculation of RE.²⁴ In contrast to the VB approach to measure aromaticity, in the block-localized wave function (BLW) approach, all basis functions and the total electrons are partitioned into few subgroups.²⁵ For example, the wave function for the benzene molecule is block localized to represent hypothetical 1,3,5-cyclohexatriene, where the double bonds resemble to those in ethylene. The energy difference between this BLW state and a real state, which is a superposition of all possible resonance states, is defined as the energy gained by the delocalization of the localized constituents (block) comparable to the ASE.²⁶ The taxonomy of energy stabilization due to cyclic electron delocalization of the π-electron includes several terms such as RE, ASE or DE with a little difference in their notion. Hence, to avoid any confusion, we choose to stick with the term ASE throughout this text. Nonetheless, both the VB and BLW approach rely on a restriction of variational space, in order to define the subunits for which the energies are calculated and compared to that of a real system. Since ASE itself is not physically observable,^26(a) direct experimental determination of ASE also becomes difficult. However, starting from the pioneering work of Pauling and Sherman,²⁷ many efforts have already been exerted to estimate ASE from experimental data. Introduction of isodesmic²⁸ and homodesmotic reactions,²⁹ has made it possible to calculate the ASE from the energies of reactants and products, obtained experimentally from calorimetric measurements.³⁰ There are two basic approaches to measure aromaticity in terms of thermochemical data.^1(b) In the first, the heat of formation of a conjugated system is compared with that of a nonconjugated model and the difference is taken as the resonance energy. However, in this method, the ASE of benzene is found to vary in the range of 22 - 64 kcal/mol depending upon the choice of suitable nonconjugated model. In another approach, the change in energy resulting from addition of hydrogen to one of the double bonds which eventually breaks the cyclic conjugation is used as the measure of ASE. The existing methodologies to measure ASE make use of the energy comparison of the hypothetical localized state and delocalized ground state of aromatic systems.^7,31 Compared to the existing methodologies, in the present work, the extra stabilization gained by the aromatic systems due to cyclic electronic delocalization is derived as a second order perturbation energy through employing an appropriate Hubbard Hamiltonian on the reference zero-order state of the aromatic system, which is constructed by localizing π-electrons on individual atoms of the ring (vide infra). This strongly localized zeroth order reference state is similar to that opted by Mo et al. in his BLW approach^25,26 or that considered by Malrieu and co-workers in the derivation of new energy based criteria Aromatic Cyclisation Energy (ACE).^19,32 They also derived the electron delocalization energy through several perturbative contributions due to delocalization of strongly localized bond molecular orbitals corresponding to a Kekule structure.^19,32 In our formalism such a model with localized π-electrons is only used to estimate inter-site charge transfer. The other parameters required to derive the energy stabilization due to π-electron delocalization is estimated on the basis of energy-minimum structures of aromatic systems. The present method avoids energy comparison between different fragments and thus gets rid of the empiricism inherent in the design of an appropriate reference state.³³ 

Benzene, being the emblem of aromatic systems, is chosen as reference system to frame a theoretical platform for measuring aromaticity in terms of ASE. A proper definition of delocalization energy requires a reference to localized picture.³² Such a model is obtained by localizing unpaired p_z electrons on each individual carbon atom of benzene (Fig. 1). In this model, the neighbouring antiparallel p_z-spins are in favourable condition to spread into one another’s non-orthogonal overlapping orbitals stabilizing the whole system. To get ASE in conjugated molecules on the basis of Pauling’s definition, Mo has used a similar model where the electronic movement is restricted within a subspace of a molecule.²⁵ The stabilization energy due to cyclic electron delocalization, defined as ASE here, adds to the ground state energy as second order perturbation contributing to the stability of aromatic system. This is a similar theoretical model as was proposed by Anderson to account for the superexchange interaction of spins.³⁴ Previously, we have also used the same model to quantify and interpret exchange coupling constant in magnetic molecules.³⁵ To conform the model constructed in Fig. 1, the system is assumed to be composed of six π-electrons, each of which is localized on the six p_z orbitals of carbon. The many electron wave function, representing the ground state of a system, can be expanded in terms of the single particle wave functions which in turn are the orthonormal spin-orbital components of a Slater determinant. This complete orthonormal set of spin-orbitals can be expressed in terms of single particle Wannier function ϕ_iσ(r − r_i), which is related to the field operator ψ_σ,^34,36

where, the f_iσ(r_i) destroys a σ spin on the spin-orbital ϕ_iσ. Since, the present work aims at estimating ASE due to delocalization of π-electrons, only the p_z atomic orbitals are used to construct the Wannier Functions. Thus the reference wave function is expanded in the basis of single particle Wannier function (WF), which can be identified as the p_z atomic orbitals on each of the carbon atoms in benzene. This approach is similar to the Huckel’s Molecular Orbital (HMO) model,³⁷ where a π–σ separation is assumed. This approximation is further extended in a different work, where the wave function of an excess electron in a system is expressed as a superposition of only the lowest unoccupied molecular orbitals of the fragments, into which the system is decomposed.³⁸ This Fragmented Molecular Orbital approach has already been successfully used to derive the valence orbital energy.

According to the Anderson’s approach, the resonating structures are connected by the hopping integrals t_ij, which carries a spin from the site i to j without spin flip. The parameter t_ij in the Hubbard Hamiltonian is obtained through the second quantization of the one-electron Hamiltonian operator (H) with the basis of Wannier functions as,³⁹ 

Operating the transfer integral t_ij twice between any two neighbouring sites to return into the initial state, and using another Hubbard parameter U, the onsite repulsion energy, Anderson obtained the stabilization of energy (ΔE) through second-order perturbation of the energy,³⁴ 

where, $f i σ †$⁠, f_iσ are the fermionic operators, creating and destroying an σ-spin at i^th site. Expanding the summation for all spin types (σ, σ′ = α, β), and using the following anticommutation rule and identities for fermionic operators,

the Eq. (3) can be transformed into,

Now, using the Jordan-Wigner transformation⁴⁰ as follows

and the anticommutation rule in Eq.(4); Eq.(5) transforms into,

Hence, through Eq.(7), the π-electron delocalization energy relates to the inter-site hopping integral (t_ij), onsite repulsion energy (U), and the up-spin and down-spin populations in the individual sites.

To get the aromatic stabilization energy according to Eq.(7), it becomes a prerequisite to organize proper computational framework for the determination of individual parameters in the equation. Prior to the computation of the parameters in Eq.(7), the geometries of benchmark aromatic systems viz., benzene, pyridine, cyclopentadienyl anion (Cp⁻), furan, pyrrole and thiophene are optimized using the PBEPBE exchange-correlation functional,⁴¹ which performs with good accuracy for wide variety of systems.⁴² However, computation of specific property in a certain class of systems is sensitive to the use of DFT functional. So, we also employ SVWN functional corresponding to local spin density approximation (LSDA). In all the cases 6-311++g(d,p) basis set has been employed. The Gaussian 09W suite of quantum chemical package⁴³ is used for geometry optimization, computation of onsite repulsion energy (U) and populations of spin orbitals (n_iσ).

Anisimov et al. have proposed to estimate the coulombic repulsion (U) between two highly correlated d-electrons as⁴⁴ 

In this work, we follow a similar approach to get the repulsion between two p_z electrons on the same site. The energy difference between (n + 1), (n − 1) and n electron states of an atom is obtained through stabilization of wave function of each electronic state. Obtained energy differences for constituent atoms of six above mentioned aromatic systems are used as U in Eq. (7). However, for heterocyclic aromatic system such as pyridine, the onsite repulsion energy cannot be same on the carbon and nitrogen p_z orbital. To solve this problem, while calculating the contribution of hopping between carbon and nitrogen, we consider the U term in Eq.(7) as an average (U_avg) of the individual repulsion energies on carbon (U_CC) and nitrogen (U_NN)

Similar problem arises in case of aromatic systems such as cyclopentadienyl anion (Cp⁻), where a lone pair participates in the π-electron delocalization. In systems as represented in Fig. 2, the π-bonding p_z orbital on the site 2 contains two electrons. Hence for delocalization, the p_z- electrons on neighbouring sites 1 and 3 have to transit through the virtual orbitals (i′) of atom 2, which will cause the value of U on site 2 to deviate from that on other sites. To consider this, the onsite repulsion energies to be used along with t₁₂ and t₂₃, are approximated as the average of individual U’s on site 1 and 2, i.e.,

In order to realize the model aromatic system with localized π-electrons in the atom-centred p_z orbitals, the approach of Marzari and co-workers are adopted and implemented through OpenMX 3.6 quantum chemistry code.⁴⁵ In their approach, the single particle WFs are assumed to be superposition of Bloch functions |ψ_n,k〉 with wave vector k and band n,⁴⁶ and can be represented as

This band (n) is considered as an isolated band, separated by a gap from the bands below and above all k points. The integration is performed over the periodic Brillouin zone (BZ), which contains all the k vectors. The effect of translating the real-space WF by a lattice vector R is taken into consideration by inclusion of the phase factor e^−ikR and V refers to the real-space primitive cell volume. In case of a discrete uniform k mesh instead of continuous k space, it is possible to have the inverse transform of the Eq. (11) as

Compared to the delocalized nature of Bloch functions, the localized Wannier functions offer a better intuitive physical picture of chemical bonding and other local correlations.⁴⁷ Marzari and Vanderbilt developed a method to localize Wannier function by minimizing the spread of density probability around localization centres⁴⁸ and the WF generated in this way is termed as maximally localized Wannier function (MLWF).⁴⁶ The advantage of using MLWFs is that they can be chosen to span a precise energy range and naturally incorporate hybridization and bonding, appropriate to their local environment (see  APPENDIX for explanation). The hopping integral between the m^th MLWF, φ_m(0), in the home unit cell and the n^th MLWF φ_n(R) in unit cell R can be expressed as,

Using the relation in Eq. (12), the Hamiltonian matrix element at a given k can be represented by,⁴⁶ 

However, by definition the Bloch states |ψ_nk〉 are gauge – invariant which induces non-uniqueness in the construction of Wannier functions from Bloch states. Hence, in order to smoothen the |ψ_nk〉 eigenstates in k– space, a unitary matrix element U_k is introduced as follows.

There are many ways to choose the unitary element. One way to do this is to use a projection technique as discussed in the  appendix. The resulting Bloch state $ψ nk ′$⁠, obtained through unitary transformation [Eq. (A2)], remains regular at all k and leads to well localized WFs.

In the present calculation, the atom-centred projections of the p_z orbital are used for constructing the π, π^* band of systems, which in fact are responsible for aromaticity. Inclusion of a range of unoccupied eigenstates allows the MLWFs to be more localized.⁴⁹ Finally, the hopping integrals between these MLWFs are calculated using linear combination of pseudo-atomic orbital (LCPAO) basis functions and pseudopotentials within generalized gradient approximation (GGA).

In the present theoretical framework, the N-electron wave function corresponding to aromatic system is described as a set of orthonormal one-electron functions. Eventually, these one-electron spin-orbitals, intrinsic to the many electron wave function, have been defined as natural orbitals by Lowdin.⁵⁰ The local 1- or 2-centre orbitals with occupancies sufficiently close to 2.00 correspond to Lewis structure and serve well to represent “true” natural orbitals for describing many electron wave function.⁵¹ To get the population of these spin orbitals, the programme NBO 3.1 of GAUSSIAN09W suit of quantum chemical package is used.⁵² In the present problem the occupancy of π-bonding orbitals is to be obtained with a provision for delocalization of π-electrons, as required to maintain aromaticity. However, the effect of delocalization appears as a weak departure from the ideal localized picture in NBO analysis and thus a better computational tool is required to properly account for the degree of delocalization of π-electrons. In this regard, we find the electron localization function (ELF)⁵ as an appropriate solution of the problem. The region of space with high value of ELF is considered to be the region with highest probability of finding electron and is called basin.⁵³ A localized domain at a given isosurface of the ELF may contain more than one such basin, which are gradually separated with increase in isovalues. The point of bifurcation is interpreted as the interaction between neighbouring basins and the isosurface value where they bifurcate is taken as a measure of delocalization.⁵⁴ In fact the ELF value has been partitioned into the σ- and π-components to constitute a new scale of aromaticity, caused by the circulation of σ- and π-electrons.⁵⁵ Hence, the ELF–π value, being a standard measure of π-electron delocalization, is multiplied with the spin orbital occupation, i.e. n_iσ and n_jσ′ (obtained through NBO) to get the actual number of delocalizing π-electron population, contributing to ASE. Since, all the systems considered in the present work exist in the singlet state, the population of a particular orbital does not alter with change in the spin part, i.e. in every occasion we get n_iσ = n_iσ′. Considering this fact and including the delocalization factor in terms of ELF–π value, the Eq.(7) takes the final form as

The ELF–π value has been obtained by employing the Multiwfn code⁵⁶ on the molecular orbitals of aromatic systems with optimized geometries.

To apply the present formalism [Eq.(16)], six different aromatic systems, benzene, pyridine, Cp⁻, furan, pyrrole and thiophene are chosen and the computed ASE values of these systems are compared with available reports.

In case of benzene, the onsite repulsion energy of a p_z– electron in carbon atom (U_CC) is calculated by applying Eq. (8) to ^•CH₃ radical, where the single spin resides on the p_z orbital. The U_CC value for each carbon atom is reported in Table I. The hopping integrals (t_ij) among six MLWFs, as shown in the Fig. 3(a), are found to be equal for any two neighbouring carbon atoms (Table I). The populations of α-spin and β-spin p_z orbitals (n_iσ), localized on each carbon atom and engaged in π–bonding, are also found to be equal for all the carbon atoms (Table I). The ELF–π value or the isovalue where the bifurcation in the π-electron cloud just begins, is found above 0.90 irrespective of functionals (Table I). This indicates almost full delocalization of the π–electrons in benzene. All these individual parameters taken together in Eq.(14), yield ASE which is quite close to the experimental value of 36 kcal/mol based on heats of combustion and hydrogenation (Table I).³³ 

In this case, the hopping integrals between MLWFs centred on carbon and nitrogen (Fig. 3(b)) vary slightly from those between neighbouring carbon pairs in the molecule (Table I). To get the onsite repulsion energy on nitrogen (U_NN) and carbon (U_CC), Eq.(8) is applied on ^•NH₂ and ^•CH₃ radicals. The average of these U_CC and U_NN values (U_avg) [Eq.(9)] is used along with t₁₆ and t₅₆. The ASE of pyridine is obtained as ∼31kcal/mol, which is within the range of reported values based on heat of combustion.³³ 

The aromaticity of cyclopentadienyl anion differs from that of benzene or pyridine in the sense, that in Cp⁻ along with the single p_z-electrons, a lone pair of electrons also participates in maintaining the cyclic electron delocalization. The value of onsite repulsion energy in the carbon atoms with singly occupied p_z orbital (labelled 1 – 4 in Fig. 3(c)) differs with respect to the carbon with lone pair of electrons (labelled 5 in Fig. 3(c)). These two types of U are obtained by employing Eq.(8) on the model systems ^•CH₃ and ^••CH₃ respectively and reported as U_CC and U_lp in Table I. The average of these two repulsion energies (U_avg) is used along with t₄₅ and t₃₅.The hopping integrals among five MLWFs (Fig. 3(c)) are found to be equal for any two neighbouring atoms (Table I).The reported values of ASE through different theoretical and experimental models are found to be self-contradictory for Cp⁻. Based on the heat of formation of benzene and that calculated by using MINDO/3 for the reaction of Cp with the anion derived from (E) -1,3-pentadiene, the ASE of Cp⁻ has been derived as 13 kcal/mol.⁵⁷ On the other hand, Bordwell et al. estimated the ASE of Cp⁻ as 24 – 27 kcal/mol on the basis of equilibrium acidity data.⁵⁸ However, the value of ASE for Cp⁻ in our approach seems to be in parity with reported values (Table I).

Computation of the onsite repulsion energy on individual atomic sites of furan follows the same strategy as opted for Cp⁻. The onsite repulsion energy on the carbon atom (U_CC) with singly occupied p_z orbital and that on the oxygen atom (U_lp) with doubly occupied p_z orbital are calculated through Eq.(8), taking ^•CH₃ and ^••OH₂ as the model systems. The average of these two values (U_avg) is used along with t₁₅ and t₄₅, i.e. the hopping integrals between the MLWFs centred on oxygen (labelled 5 in Fig. 3(d)) and neighbouring carbon atoms (labelled 1 and 4 in Fig. 3(d)). Though the estimated value of ASE in GGA functional is quite close to the reported values,^14,27,33,59 this agreement is not found in case of LSDA functional (Table I).

The onsite repulsion energies on the carbon atoms (U_CC) with singly occupied p_z – orbital and nitrogen atom (U_lp) with doubly occupied p_z orbital are calculated by employing Eq.(8) on model systems ^•CH₃ and ^••NH₂ respectively. The average of these two values (U_avg) is used along with the hopping integrals (t₁₅ and t₄₅) in between nitrogen and its neighbouring carbons. The ASE of pyrrole, calculated by Pauling-Sherman method⁶⁰ using Pauling’s table of bond-energies and thermal data from Landolt-Bornstein has been reported as 23 kcal/mol.⁶¹ Whereas, the homodesmotic reaction scheme of Cyranski et al. produces 21 kcal/mol ASE for Pyrrole.¹⁴ On the basis of heat of combustion, the ASE of Pyrrole has been reported to be 14 – 31 kcal/mol.³³ Hence, the ASE 24 kcal/mol in GGA functional and 15 kcal/mol in LSDA functional obtained through Eq.(14) lie within the range of reported values (Table I).

The onsite repulsion energies on the carbon and sulphur atoms are evaluated through implementation of Eq.(8) on the optimized model systems ^•CH₃ and ^••SH₂. The values are reported as U_CC and U_lp respectively in Table I and the average of these two values (U_avg) is used along with t₁₅ and t₄₅, the hopping integrals between the MLWFs centred on sulphur (labelled 5 in Fig. 3(f)) and its neighbouring carbon atoms. Compared to furan and pyrrole, the longest bond-length between sulphur and carbon in thiophene causes the hopping integral between these two to be the lowest (Table I). However, for the lower electronegativity of sulphur in thiophene compared to oxygen in furan and nitrogen in pyrrole, the π-electron can be maximally delocalized in thiophene. This fact is attested through the highest ELF-π value of 0.81 in thiophene in GGA functional, though the LSDA functional gives a much underestimated value of ELF-π. In accordance with the homodesmotic reaction scheme proposed by Cyranski et al., the ASE of thiophene has been reported as 19 kcal/mol,¹⁴ which is close to the present estimation of 22 kcal/mol (Table I). However, the ASE for thiophene is reported as 31 kcal/mol through Pauling–Sherman method,⁶⁰ 24 – 31 kcal/mol on the basis of heat of combustion and 29 kcal/mol based on heat of hydrogenation.³³ 

This work exerts an effort to quantify aromaticity, which has long been a challenging task in the domain of Chemistry. Till date, a plethora of aromaticity indices have been introduced as an indirect measure of aromaticity due to its multidimensional nature. In spite of the inherent diversity, the aromaticity indices are found to stem from a single phenomenon, which is the cyclic delocalisation of electrons. Hence, in the present work, a proper quantum chemical protocol is figured out to quantify aromaticity in terms of the energy stabilization due to cyclic delocalization of π-electrons, termed as the Aromatic Stabilization Energy (ASE). In the proposed phenomenological model, the ASE is equated with the second order correction to the reference state energy resulting from atom to atom charge transfer. The expression of the second order perturbation energy due to atom to atom charge transfer is obtained on the basis of Anderson’s model of superexchange. Subsequent modification of the expression gives the final form of ASE in terms of the inter-site charge transfer integral, onsite repulsion energy and population of spin orbitals. To estimate the inter-site hopping integral, a reference state of the aromatic systems is constructed where the neighbouring atoms of the ring are considered to have antiparallel oriented π-electrons, localized at their respective p_z orbitals. Such a reference state is achieved in the form of atom-centric Wannier functions with the character of p_z orbital. Thus, the compartmentalization of the total molecular space into atom-centric maximally localized Wannier functions enables estimation of inter-site charge transfer integrals. On the other hand, the onsite repulsion energy (U) at each atom of the ring is evaluated as the energy difference between the neutral and ionic states of a molecular fragment with similar electronic environment to that atom. In order to get the exact occupancy of an atom-centred p_z spin orbital, half the population of π-bonding orbital is multiplied with the ELF-π value, which in essence is the measure of electron delocalization. Taking the degree of electron delocalization into account, it becomes possible to derive the exact population of the spin-orbital maintaining the ring current and thus contributing to the ASE. Estimation of the hopping integral, onsite repulsion energy and the population of the spin-orbital through different computational schemes and using these parameters in Eq. (16) ultimately give the value of ASE, a well established direct measure of aromaticity. To verify the general applicability, the devised formalism is employed on six different benchmark aromatic systems. The values of ASE obtained in all the cases, particularly in GGA functional, are found to be concordant with available experimental and/or theoretical estimates, which in consequence assure the validity of present formalism. However, to establish the general nature of the develop formalism, it becomes urgent to apply the present approach on polyacenes and novel all metal aromatic systems.⁶² 

This work is dedicated to Professor S. N. Datta of Indian Institute of Technology, Bombay on the occasion of his sixty fifth birthday. Financial support from the Council of Scientific and Industrial Research (CSIR), India is gratefully acknowledged. TG thanks CSIR for the senior research fellowship.

To construct a smooth gauge in k-space and corresponding set of well-localized WFs, a projection technique is adopted, where one starts from a set of J localized trial orbitals, ϕ_n(R) corresponding to rough guess for the WFs in the home unit cell. Within the continuous k formulation, these ϕ_n(R) s are projected onto the Bloch function at wave vector k to obtain

In practice, such projection is achieved by computing a matrix of inner products (A_k)_mn = 〈ψ_mk|ϕ_n〉 followed by the overlap matrix $( S k ) mn = 〈 ϕ mk | ϕ nk 〉 V = ( A k ∗ A k ) mn$⁠. These are next used to construct the Löwdin-orthonormalized Bloch-like states

where $| ψ nk ′ 〉$ has now a smoother gauge in k-space. When these are put back into Eq. (11) in place of |ψ_nk〉, it results in a set of well-localized WFs, which naturally incorporates hybridization and bonding appropriate to their local environment.