On the origin of bonding in metals: lithium as a case study

The bonding in metals is analysed within the framework of the PATMOS (Perturbed AToms in MOlecules and Solids) model. The electronic binding energy per atom is written as a sum of a distortion energy of the atoms and the partitioned interaction energy comprising Coulombic, exchange and correlation terms. The adopted physical model of the infinite system, is spherical embedding of the atoms of the reference unit cell. Correlation energies are calculated by second-order Møller-Plesset and second-order Epstein-Nesbet perturbation theory. The binding energy of lithium solid is calculated for 16 nearest neighbour distances from 4.0 to 10.0 Bohr. Electron correlation is of paramount importance for the binding energy. The calculated cohesive energy is Hartree comparing with experimental value 0.0599 Hartree. The bonding picture is characterised by slightly expanded atomic orbitals. GRAPHICAL ABSTRACT


Introduction
Metallic bonding is apparently very different from bonding in polyatomic molecules. The latter is localised and directional. For a huge class of molecules, it can be described by localised electron-pair bonds [1]. According to Ruedenberg and coworkers [2], it is the quantum mechanical interference of wave functions of fragments that leads to the kinetic energy decrease in the system -CONTACT Inge Røeggen inge.roeggen@uit.no Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, UiT The Arctic University of Norway, Tromsø N-9037, Norway which becomes the driving force behind chemical bonding. Most of studies have agreed upon the Ruedenberg's framework though fine details may diverge for different theoretical methods [3]. This picture of electronpair bonding is not easily transferred to a metal. For an alkali metal with a bcc-structure, there are eight nearest-neighbour atoms for a particular atom, but only one valence electron for each atom. Hence, a system of electron-pair bonds for the metal is indeed difficult to imagine. To cope with this enigma it is customary to suppose that a metal is characterised by a delocalised 'sea' of electrons and that this 'sea' of negative charge holds the atoms together. The corresponding mathematical model is Hartree-Fock theory with delocalised orbitals, i.e. canonical Hartree-Fock orbitals. However, within the one-determinant approximation, the wave function is invariant with respect to a unitary transformation of the orbitals. A delocalised picture can be transferred to a localised picture. Accordingly, the original problem of explaining bonding is as puzzling as before. The scientific community has then more or less disregarded this particular problem and instead focussed on what can be obtained by molecular orbital theory with delocalised orbitals. An extensive body of scientific work demonstrates the success of this approach.
One of the difficulties of explaining bonding in metals can be traced back to a simplistic interpretation of the electron-pair bond concept. It is often interpreted as if electrons in complexes prefer to stick together in pairs. However, electrons repel each other, and left to themselves they separate. As for the ground state of a system, electrons try to come as close as possible to the nuclei in accordance with the Pauli exclusion principle. For localised orbitals we might interprete the Pauli principle as stating that two electrons with different spins can occupy the same part of the physical space. In order to obtain the lowest possible total energy for the system, the electrons arrange themselves in such a way that they are close to a nuclues or a pair of nuclei. In electron-pair bonds an electron of an atom is shifted towards the nucleus of a neighbouring atom, and vice versa. This particular feature of the electron-pair bonding is demonstrated in a work by Røeggen and Gao [4]. By moving away from the paring concept but keeping its physical content of the bonding, i.e. shifting of electrons in direction of neightboring atoms, bonding in metals might be more easily understood.
Another key to understand bonding in metals is related to the electron density. The density in metals is not very different from the sum of the atomic densities of isolated atoms put in the positions of the nuclei of the metal. Hence, the atomic wave-functions are only modestly distorted in forming the metal. This fact has an important corollary. If one uses restricted Hartree-Fock in the study of metals, which implies doubly occupied orbitals, this feature of the metallic bond cannot be disclosed. Therefore, one should preferably use unrestricted Hartree-Fock as a basic approximation in describing metals since it is essential to have one spatial orbital for each valence electron.
The purpose of this work is twofold: first, to construct a workable computational model for ab initio studies of the electronic structure of atoms in three-dimensional lattices. Second, to understand the origin of bonding in metals. As for the first purpose, the PATMOS model [4,5] (Perturbed AToms in MOlecules and Solids) is a convenient starting point. PATMOS calculations on finite chains of metallic atoms strongly suggest that the electronic structure of the atoms of chains can be characterised by localised orbitals, i.e. distorted atoms. Hence, a physical model with spherical embedding around each atom of the reference unit cell, is appropriate. The spherical embedding guarantees correct symmetry. In this work we shall then demonstrate that ab initio calculations on metallic atoms in three-dimensional lattices are feasible.
The structure of the paper is as follows: The second section is devoted to a short description of the PATMOS model adapted to periodic systems. The third section focuses on a physical model of the solid using spherical embedding. The fourth section gives details on the computational aspects of the model. The fifth section is devoted to calculations on lithium atoms in a threedimensional lattice. In the last section there is a discussion of the different sources of errors in this work.

The PATMOS model
To manage the electron correlation problem for large molecules, it is now well established that within the wave function framework local correlation models are required. The local models can be traced back to the fundamental works by Pulay and coworkers [6][7][8][9]. To avoid a large virtual orbital space, Pulay [6] proposed to remove the component of the localised occupied orbitals from the atomic basis. He denoted the modified basis projected atomic orbitals (PAO). In calculating the inter-atomic correlation energy one could then use the union of the PAO's for the atoms involved instead of the full virtual orbital space for the complex. Hence, a drastic reduction of the computation time could be achieved.
To date various local correlation models have been developed and applied for periodic systems, for instance, the local second-order Møller-Plesset theory (MP2) [10], divide-expand-consolidate MP2 [11] and clusterin-molecule local correlation approach [12]. Another way of adapting local correlation methods to metals is to use the energy incremental scheme. The orbital energy incremental scheme was introduced by Nesbet [13], and refined by Stoll and coworkers [14][15][16], Røeggen [17,18], Bytautas and Ruedenberg [19], Voloshina and Paulus [20], and more recently by Paulus and Stoll [21]. The PATMOS model is formulated within the energy incremental scheme. In principle the correlation energy in PATMOS model can be calculated by any size extensive correlation model. The detailed description of the PAT-MOS model is given in our previous works [4,5]. In this work we sketch only the essential elements: the energy partitioning and the PATMOS basis set approach.

Energy partitioning
The total electronic energy of a complex can within the PATMOS framework, be written as In Equation (1) where An important term for a periodic system is the sum of effective atomic energies for the atoms in the reference cell: In Equation (4) N uc atoms is the number of atoms in the unit cell. As we increase the number of atoms in the physical model, E PATMOS uc should converge to the corresponding value for the infinite system.
In this work we are interested in the binding energy per atom: where E iso A is the energy of the isolated atom. Then it follows from Equation (3) where dist The term dist A , the distortion energy, represents the change in the energy of atom A in the complex due to the presence of the surrounding atoms. The interpretation of the terms in Equations (8) to (10) should be evident.
The virial theorem yields additional insight into the character of electronic binding energy. Let T denote the kinetic energy operator: where n elec comp is the number of electrons in the considered complex. Then according to the virial theorem: where E kin iso,exact is the kinetic energy of the complex comprising the isolated atoms or subsystems, calculated by the exact wave function. A similar reltion is valid for the bonded complex with the equilibrium geometry: Hence, we have for the binding energy: Accordingly, if the system is bounded, then E kin comp,exact > E kin iso,exact .
As a consequence, in the energy partitioning scheme the change in the kinetic energy is a 'response' effect, but it is hidden in the positive distortion energy.

The PATMOS basis set procedure
In the recent work by Røeggen and Gao [5] a basis function (BF) region is defined as a unit cell and a certain number of nearest neighbour unit cells. However, since we in this work advocate spherical embedding, the BFregion of an atom comprizes the atom in question and all its partner atoms in a sphere centred on the nucleus of the atom. The number of atoms of the BF-region is denoted N A BF . See illustraction in Figure 1. Two different atomcentred basis sets are associated with each nucleus, a large one, {χ A,lb μ ; μ = 1, . . . , m lb A }, and a small one, {χ A,sb μ ; μ = 1, . . . , m sb A }. The basis set for an atom in the unit cell is then where A BF is the set of atoms in the BF-region (excluding atom A).
The spatial part of a spin orbital (α-or β-type): The orbitals of the UHF wave functions are subjected to orthogonality constraints: The constraints, Equations (18) and (19), require special attention in the optimisation procedure since the orbitals of different atoms are expressed in different basis sets.

Orbital localisation measures
To describe orbital localisation, we use charge centroids and charge ellipsoids [22] in this work. As illustrated in our previous work [4], it is a very compact and visual way of looking at orbital localisation. The charge centroid is defined as where ψ is any spatial orbital. The extension or spread of an orbital can be described by the second-order variance matrix where x C r is the rth component of the charge centroid vector r C . Diagonalization of the second-order variance matrix M rs yields the charge ellipsoid.
The eigenvalues {a 1 , a 2 , a 3 } of the matrix M rs correspond to the squares of the half-axes of the ellipsoid. The standard deviations in three orthogonal directions, i.e. the directions of the half-axes, are therefore given by which can be used as a measure of the extension of the orbital with respect to the charge centroid position. We may also use the volume of the ellipsoid V = 4 3 π l 1 l 2 l 3 as a single number for the extension.
The half-axes also define a lower bound for the kinetic energy associated with a given orbital ψ: (23) This inequality can be derived from the Heisenberg uncertainty principle [4] and expresses neatly that the kinetic energy increases when an orbital becomes more compact.

The physical model for a periodic system
In the physical model for constructing the orbitals of the UHF function we include as a constraint that all atoms of the same type should be described by identical orbitals, i.e. translational symmetry should be satisfied. Hence, there are no surface effects in the physical model. In Figure 1 we consider a one-dimensional system with two different atoms in the unit cell. We look at two different types of embedding. The first one, (a), is an embedding with a fixed number of unit cells. In Figure 1, embedding (a) consists of the reference cell and the two first nearest neighbouring cells (1NN). The second embedding is a spherical one including all atoms in sphere surrounding an atom of the reference cell. In Figure 1, the radius of the sphere is slightly larger than the unit cell distance.
When we optimise the UHF orbitals, we consider an effective Hamiltonian for each atom in the reference cell. The Hamiltonian includes only the atoms of the BFregion. We notice that for embedding (a), the atom A of the reference cell has two atoms on the left and three atoms on the right, and vice versa for atom B. This lack of symmetry creates an artificial shift of the charge density. For a one-dimensional system this shift can be reduced by including a large number of unit cells in the BF-region. However, for atoms in three-dimensional lattices, such an extension is not feasible. On the other hand, if we look at spherical embedding, Figure 1(b), we have perfect symmetry for any size of the embedding.

Computational aspects
The computational feasibility of a model is of paramount importance when we consider large complexes such as a metal. In the first subsection we shall give details on the procedure for obtaining the spin orbitals of the UHF wave function. The second and the third subsections are devoted to the electron correlation problem restricted to intra-atomic and diatomic terms. In particular the diatomic terms need careful considerations. We have to compute of the order 1000 diatomic terms for obtaining sufficiently converged results. Hence, a simplistic, but still fairly accurate, procedure has to be constructed.

Determination of the UHF orbitals
The starting point for the UHF procedure is an effective Hamiltonian for an atom A of the reference unit cell: where the effective one-electron Hamiltonian is given by N A and N B are the number of electrons of atoms A and B, respectively. J B j and K B j are the Coulomb and exchange operators generated by the spin orbital ψ B j of atom B, and h(r i ) is the one-electron operator for the complex: atom In Equation (26), the symbols have their conventional meaning. The spin orbitals of the atom A are orthogonal to all spin orbitals of the atoms in the set A BF , i.e.
In order to comply with the orthogonality constraints of Equation (27), and the translational symmetry of the orbitals, the optimisation procedure consists of several iterative steps. Preliminary calculations on one-dimensional systems of pure metals demonstrate that the localised orbitals are slightly perturbed atomic orbitals. The PATMOS basis set procedure with a large and a small basis sets assigned to each atom, is particularly well-suited to deal with this feature. The large basis set can easily account for the small expansion of the atomic density while the small basis set for the partner atoms in A BF takes care of the orthogonality constraints. However, we can introduce a further simplification. Let P B,α occ denote the projection operator associated with the α-type spin orbitals of an atom B in the BF-region of atom A. The dual space of atom A, Equation (16), is modified in the follwing way: The next step is to solve the Hartree-Fock equations for atom A: and The orbitals are expressed in terms of the modified large basis sets, i.e. and A fixed-point iteration procedure does not necessarily converge. Hence, we introduce a UHF functional where The orbitals {ψ A i (λ)} have to be properly symmetrised. The procedure is described in Appendix. We denote the symmetrised orbitals {ψ then we continue to the next iteration. Otherwise, we choose new values of λ until we can perform a parabolic fit. Two or three values of λ are in most cases sufficient.
To obtain a converged result of accuracy 10 −7 Hartree, typically 15-20 iterations are required.

Correlation energy
Typically, intra-atomic MP2 energy yields 80-90% of the corresponding full configuration interation (FCI) energy. The intra-atomic MP2 energy for lithium varies only slightly for different values of the lattice parameter, i.e. less than 1%. Hence, MP2 should be sufficiently accurate for our purpose. The inter-atomic correlation terms are more demanding. In this work we are using both MP2 and secondorder Epstein-Nesbet (EN2) theory. MP2 underestimates the inter-atomic correlation enery while EN2 usually overestimates it.

Intra-atomic correlation energy
where the two sums are restricted to respectively occupied and virtual spin orbitals. By distinguishing between αand β-type spin orbitals, the MP2 energy can be written as a sum of three different terms:

Inter-atomic correlation energy
For atom A we calculate the diatomic correlation terms for all atoms with a sphere of radius r sph model surrounding the nucleus of A. The radius of the sphere is determinde by a selected convergence criterion. For a pair of atoms (A, B) we define inter-atomic embedding A inter and B inter . The embedding A inter includes all atoms within a sphere of radius r sph inter surrounding the nucleus of atom A (but excluding atom A). Similarly for atom B. The embedding for the pair (A, B) is then (see Figure 2) The effective Hamiltonian for the pair (A, B) is where In Equation (41), the symbols have their standard meaning. The Coulomb and exchange operators in Equation (41) require special attention. For each pair of atoms (A, B), the proper calculation of two-electron integrals requires three sets of dual basis sets as defined in Equation (16): A dual , B dual and C dual . Since the number of atoms in AB inter is of the order of the sum of atoms in the sets A BF and B BF , the computational cost of this proper approach is prohibitively high. Our simplistic procedure is based on the fact that the charge ellipsoids [4] associated with the occupied orbitals are spherical symmetric. Hence, we project the orbitals {ψ C i ; 1 ≤ i ≤ N C } onto the small basis set {χ C;sb μ ; 1 ≤ μ ≤ m sb C }:  (27)]. If B / ∈ A BF , then the orthogonality condition is not necessarily fulfilled. Consequently, we perform a symmetric orthonormalization of the occupied orbitals. This is followed by a Gram-Schmidt orthonormalization of occupied plus virtual orbitals. The orbitals are appropriately divided in three separate groupsfor α-type orbitals, {ψ A;occ;α By using spin orbitals the diatomic MP2 energy can be written as and h AB eff (r) is defined in Equation (41). The Coulomb and exchange operators in Equation (50) have their conventional meaning.
By introducing the following set of two-electron functions In actual calculations the expression for EN2, Equation (53), is slightly modified. The spin orbitals ψ A i and ψ B j are determinde by using a spherical embedding. As is evident from Figure 2

The bcc structure of lithium
In this section we present a series of calculation on solid lithium in a body centred cubic (bcc) structure (unit cell parameter or lattice parameter a). The lithium atoms situated at the corners of cubes have the electronic configuration (1sα1sβ2sα), for short α-atoms, and the atoms in the centre of cubes (1sα1sβ2sβ), for short β-atoms. The unit cell (uc) includes one α-atom and one β-atom.
The basis sets adopted are uncontracted Gaussian type functions (GTF). A small basis set (sb) = (11s1p), and two large basis sets (lb) = (15s7p2d1f ) and (lb) = (19s8p7d5f 2g1h). Our integral code requires family type basis sets. Further, the exponents are all drawn from the same set of universal type exponents, i.e. {η k = αβ k−1 , 1 ≤ k ≤ k max }. The parameters defining the basis sets are given in Table 1.
In this work two-electron integrals are represented by Cholesky vectors. The idea of a Cholesky decomposition of the two-electron integral matrix was first suggested by Beebe and Linderberg [23] and later adopted by several research groups [24][25][26]. In this work, an integral threshold of 10 −7 Hartree is used for the calculation of intra-atomic terms, and 10 −5 Hartree for inter-atomic terms. The latter approximation affects the total energy per atom by less than 10 −6 Hartree.
The BF-region is defined by a sphere of radius, r sph BF = a uc . A sphere with the cube corner (0, 0, 0) and radius r sph BF encloses six α-type atoms and eight β-type atoms. The first nearest-neighbour distance r 1NN is the distance between the α-atom in position (0, 0, 0) and the β-atom in position ( a uc 2 , a uc 2 , a uc 2 ). The UHF and MP2 energy for isolate lithium is given in Table 2. The half-axes of the charge ellipsoid of the 2s orbital are also included.
The typical pattern of convergence of the iterative procedure for obtaining the perturbed atomic orbitals, is shown in Table 3. Only atoms of the BF-region are included in the determination of the orbitals. The orbitals of the isolated atoms, properly symmetrised for the BFregion, are used as starting orbitals. The effective UHF energy for an atom is converged after 14 iterative cycles. For all cases considered in this work, convergence is obtained after 15-20 iterations. In order to calculate the interaction between an atom A in the reference unit cell and all atoms within a sphere of radius r sph model surrounding the nucleus of atom A, we have partitioned this sphere into an inner sphere of radius r = a uc , and concentric shells of thickness a uc . The results are presented in Table 4 for a nearest-neighbour distance r 1NN = 6.0 Bohr. We notice the rapid convergence with respect to shell distance. In this work shell 4 is the last shell included. The correlation terms are not fully converged, but the contribution from neglected shells are assumed to be small compared with the total error of the EN2 correlation energy.
The main results of this work is presented in Table 5 and Figure 3. Binding energy per atom as a Table 1. Lowest and highest exponents of the GTF-family basis sets used in this work (η k = αβ k−1 , 1 ≤ k ≤ k max ).  function of the first nearest-neighbour distance, r 1NN , is given for UHF, UHF+MP2 and UHF+MP2/EN2. The experimental cohesive energy [27] for lithium is 0.059903 Hartree (cohesive energy is defined as positive). The nearest-neighbour distance r 1NN = 5.744214 Bohr corresponds to experimental crystal structure [28]. We notice that UHF yields only a shallow minimum of −0.0047 Hartree at a distance of 7.03 Bohr (parabolic fit). The UHF+MP2 gives a binding energy of −0.028011 and an equilibrium distance of 5.997 Bohr (parabolic fit). The PATMOS-MP2 result for the binding energy is less than 50% of the experimental binding energy (correcting for zero-point energy). The PATMOS-EN2 binding energy is roughly 5% off the experimental value. However, the EN2 inter-atomic correlation energies are subjected to unsystematic errors. As is evident from Figure 3, the value for r 1NN = 5.4 Bohr, i.e. −0.062566 Hartree, overestimates the binding energy. If we disregard this value and perform a parabolic fit based on the distances {5.0, 5.2, 5.6} Bohr, we obtain a binding energy of −0.05817 Hartree at a distance of r 1NN = 5.30 Bohr. We shall return to this problem in the error discussion, Section 6. The magnitude of the difference between the estimated value of the binding energy and the calculated value for r 1NN = 5.4 Bohr, is used as a measure of the unsystematic error of the EN2 correlation energies. It is roughly 8% of the EN2 correlation energy.  The character of the binding of solid lithium is disclosed in Table 6. To easily group the relative contribution of the components, we have used the magnitude of the Coulomb interaction as energy unit. The partitioning of the binding energy demonstrates very clearly that the correlation energy is of paramount importance. Its relative contribution increases with increasing value of r 1NN . There is one repulsive term: the distortion energy, i.e. the charge in intra-atomic energy due to the interaction with the surrounding atoms. The magnitude of the 'semiclassical' Coulomb energy, E Coul Li;inter , is less than the distortion energy for all distances considered in this work. As a consequence, the Hartree model, i.e. the UHF energy minus the exchange energy (UHF orbitals are used for the Hartree model) yield no binding. This feature is demonstrated in Table 7. Measured by the half-axes of the charge ellipsoid of the 2s orbital, we notice that there is only a small extension of the electron density at the experimental equilibrium distance r 1NN = 5.744214 Bohr.
With exception of some round-off errors less than 10 −6 Hartree, we have identical results for the αand βatoms of the reference unit cell. As for the inter-atomic interactions, those terms are only calculated for the α-type atom.
To summarise this section: binding of lithium in a bcc structure can be described by localised orbitals. The perturbed atomic orbitals are only slightly distorted. Furthermore, the correlation energy is extremely important in order to obtain an accurate binding enregy.

Discussion
There are mainly three different sources of error in this work: basis set, correlation energy and model errors. We consider first the basis set errors. At the distance r 1NN = 6.0 Bohr, there is a small extension of the atomic electron   density, see Table 7. The half-axes of the charge ellipsoid of the 2s orbital are respectively 2.504 Bohr and 2.498 Bohr for the two large basis sets (15s7p2d1f ) and (19s8p7d5f 2g1h). As demonstrated in Table 8, there are only small differences of the UHF energy components, roughly 1% or less. For the MP2 inter-atomic correlation energy the change is around 3%. Due to the unsystematic error of the EN2 component, the basis set effect is hidden in the unsystematic error. If we consider the difference between the two EN2 energies in Table 8 as a measure of the unsystematic error in the EN2 calculations, the error is roughly 7% of calculated correlation energy. This estimate is consistent with our previous one obtained in Section 5, i.e. 8%. In Table 9, the change in intraatomic MP2 energy is shown as a function of the r 1NN distance. The change is less than 1% for all distances considered. For the smaller large basis the change in intraatomic MP2 energy due to the surrounding atoms, is 0.000060 Hartree to be compared with 0.000056 Hartree for the basis (19s8p7d5f 2g1h). Hence, the basis set error with respect to MP2 intra-atomic correlaton energy is negligible.
There are two types of model errors in this work related to the spheres with distances r sph model and r sph inter , see Figure 2. The error of neglecting atoms outside the model sphere, i.e. the larger one, is well accounted for in this work. However, the interpair spheres with radius r sph inter are the main sources of errors. In principle it can be reduced by increasing the radius r sph inter , say from a uc to 2a uc . But then the number of atoms in interpair spheres increases from 14 to 64, and the number of basis functions increases by almost a factor 5. Hence, at present this is not computationally feasible.
Our best estimate of the binding energy derived from a parabolic fit using the nearest-neighbour distances {5.0, 5.2, 5.6} Bohr: −0.05817 Hartree. The binding energy per atom as a function of r 1NN , i.e. E bind Li (r 1NN ), is describing a collective contraction/expansion of the crystal. In approximating the nuclear motion of an atom, we use a harmonic potential. A straightforward calculation yields a zero-point vibratonal energy per atom equal to 0.00076 Hartree. With a 7% unsystematic error of EN2 correlation energy we arrive at the following estimate of the cohesive energy 0.0571 ± 0.004 Hartree, comparing with experimental value [27] 0.0599 Hartree.

Concluding remarks
Our main objective of this work has been achieved: to construct a computationally feasible model for describing the electronic structure of pure metals using localised orbitals. Our calculated cohesive energy for lithium is in fair agreement with the experimental result.
A second point to emphasise is the importance of the electron correlation energy. The UHF model yields only a shallow minimum of the potential energy surface. Large basis sets are necessary conditions for obtaining reliable binding energies. A strong merit of the PATMOS model is the dual basis set approach. Large basis sets of the same quality as for molecules can be adopted. Furthermore, since we are using fixed atomic domains for the orbital spaces (slightly modified by orthogonality constraints), and we are calculating only diatomic corrections, linear dependencies can be avoided by just checking for linear dependency in calculations on the corresponding diatomic molecules.
A challenge for future work is to extend the model such that the unsystematic errors of the EN2 correlation energies can be reduced. Further, to improve the calculated correlation energy, we shall consider a more accurate model than EN2 for the nearest-neighbour atom pairs.
We are now in the beginning of a research programme devoted to the study of chemical binding using the PATMOS framework. As for solids, this implies timeconsuming calculations. Some preliminary calculations on solid beryllium seem to indicate that bonding in this case might not be correlation driven. The results will be reported in due time.
Finally, if our conjecture concerning localised orbitals and bonding in metals, can be confirmed, then there is nice similarity between molecules and metals: bonding for both types of systems is most appropriately described by localised orbitals, but for spectroscopy canonical orbitals is the proper choice.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by the Research Council of Norway (Norges forskningsråd) through its Centres of Excellence scheme [project number 262695]. This work has received support from the Norwegian Supercomputer Program (NOTUR) through a grant of computer time [grant number NN4654K].

Appendix. Symmetric orthonormalization of spin orbitals
Let A denote an atom in the reference unit cell, and A p an atom of the same type in the BF-region of atom A.
The iterative cycle is then repeated. Typically, 15-20 cycles are required in order to have a converged result. If the procedure does not converge within a fixed number of cycles, the parameter λ in Equation (35), i.e. the step length, is reduced.