text:
# Hartree-Fock-LAPW method: using the full potential treatment for exchange

## I Introduction

In contrast to successful applications of the Hartree-Fock (HF) method to molecules, which have become a routine procedure, Hartree-Fock calculations of electron band structure of solids are relatively rare. The electronic band structure calculations mainly rely on the density functional approach  with the local density approximation (LDA) or other non-local exchange approximations  which try to mimic the Fock exchange for conduction electrons. As known, LD approximation is a product of the Hartree-Fock theory of the free electron gas and as such it inherits its shortcomings.  The LD approximation being universal and simple from one side, is structureless and leads to an unphysical self-interaction,  from the other. In addition, LDA is an uncontrolled approximation, and there is no straightforward path to further improvements in its accuracy.  Therefore, it is highly desirable to combine the Hartree-Fock approach with methods based on plane wave basis sets which traditionally play an important role in the field of electron band structure calculations. Although a variety of computational procedures on periodic HF method has been reported in the literature,  all of them use Gaussian-type orbitals (GTO). In this paper we present a merger of the Hartree-Fock approach with the linear augmented plane wave  (LAPW) method, which constitutes a new Hartree-Fock-LAPW procedure for electron band structure calculations. The consideration of the LAPW basis is important since the method has been extensively tested and it is traditionally associated with electronic band structure calculations.

At the center of the Hartree-Fock-Roothaan approach there are calculations of the direct Coulomb and the exchange matrix elements  and in the present work we consider thoroughly all cases, including the delocalized and localized states. The advantage is that all expressions for the matrix elements can be found in one source (i.e. this paper), but on the other hand, this has led to numerous formulas and their extensive derivations. While GTOs are popular in calculations of molecular systems and have the advantage of giving easily integrable polycentric functions, we lose this benefit working with the LAPW basis set. There, we have to rely on a different strategy of calculations. Our choice here is the full potential LAPW (FP-LAPW) treatment,  which is based on the procedure of finding the general solution for a periodic Poisson’s equation. For Poisson’s equation we use a new original technique (Sec. II). It employs the Ewald method  for the potential of the monopole terms and an exact expansion in Fourier series in the interstitial region, with a multipole series expansion of the potential inside the spheres.

A principal new element of our work is to use the strategy of the FP-LAPW method for the calculation of exchange matrix elements between itinerant electrons. Indeed, an exchange matrix element $`ab|V^{Coul}|ba`$ ($`a`$ and $`b`$ here represent electronic states with the same spin component) can be estimated as a Coulomb self-interaction energy of the charge density distribution $`\psi _b^{}(\stackrel{}{R})\psi _a(\stackrel{}{R})`$, which appears as a result of the electron transition $`ab`$. The procedure implies (i) an introduction of an auxiliary potential of the charge density $`\psi _b^{}(\stackrel{}{R})\psi _a(\stackrel{}{R})`$ and (ii) the integration of the potential with the charge density. Such treatment is rigorous and ensures that all contributions including the long range ones (like the exchange Ewald sums) are carefully considered. It is worth to notice that a simple truncation of the exchange series of other methods excludes such long range effects.

Here we also consider crystal electric field (CEF) effects for the closed shell core electrons. Usually in the full potential LAPW treatment the core electrons are taken in the spherical approximation and the CEF effects are ignored.  While the crystal field splittings are small it may be still necessary to include these effects for precise calculations of solids where one compares the stability of different crystalline phases where the energy difference is of the order of 100 K. The present consideration is a generalization of previous works on CEF effects.

For clearness and simplicity in the following we consider only one atom per unit cell. However, the method can be easily generalized for the case of a few atoms. In this paper we consider the core states being completely confined inside the spheres. Various improvements of the method in this respect are also possible,  but they are not of fundamental importance and will not be concerned here.

The paper is organized as follows. We start with the general solution of periodic Poisson’s equation, Sec. II. In Sec. III we illustrate how this technique works for the case of a face centered cubic (fcc) crystal. We compare the convergence and accuracy of different approaches. Then we discuss the Hartree-Fock-Roothaan method for solids with extended states, Sec. IV, and calculate the direct Coulomb matrix elements for LAPW basis functions. In Sec. V we present calculations of the matrix elements of exchange, while Sec. VI is reserved for our conclusions.

## II General solution of periodic Poisson’s equation

The general solution of Poisson’s equation $`V^{Coul}`$ for a periodic charge density was demanded by necessity to improve the “muffin-tin” (MT) approximation of LAPW method.  LAPW electron band structure calculation procedure with the general potential $`V^{tot}(\stackrel{}{R})=V^{Coul}(\stackrel{}{R})+V^{exc}(\stackrel{}{R})`$ is referred to as the full potential LAPW (FP-LAPW) method  since in principle no shape approximation (except LDA) for the potential is assumed. First full potential calculations have been reported in Refs. , and since then they are often used to investigate the electron band structure of solids. These FPLAPW calculations are based on the general solution of Poisson’s equation given by Weinert in Ref. . The basic idea which has been proposed already in Ref.  includes two steps: 1) to obtain the potential in the interstices and then 2) to solve the boundary value problem inside a sphere. The point 1) in Ref.  represents an alternative to the Ewald method  where the true charge inside the spheres is replaced with a pseudo-charge-density of the same multipole moments. In the present paper we propose a new technique which uses the Ewald expansion for the monopoles ($`l=0`$) and the exact Fourier expansions for the higher multipoles ($`l>0`$) in the interstitial region. The procedure is straightforward and besides a truncation of the multipole series, it contains only one cut-off parameter for the Ewald expansion. The strategy consists of the same two steps as outlined above.

### A Dual and equivalent representations of charge density

We start with the description of charge density in a crystal. For a periodic crystalline structure we introduce $`\stackrel{}{n}`$ to label the sites and the direct lattice vector $`\stackrel{}{X}(\stackrel{}{n})`$ which specifies the centers of atoms (spheres). The position vector $`\stackrel{}{R}`$ inside a sphere $`S(\stackrel{}{n})`$ is given by

$$\stackrel{}{R}(\stackrel{}{n})=\stackrel{}{X}(\stackrel{}{n})+\stackrel{}{r}(\stackrel{}{n}).$$
(1)
In spherical coordinates $`\stackrel{}{r}(\stackrel{}{n})=(r(\stackrel{}{n}),\widehat{r}(\stackrel{}{n}))`$, where $`\widehat{r}\mathrm{\Omega }=(\theta ,\varphi )`$. For a Bravais lattice of $`\stackrel{}{X}(\stackrel{}{n})`$ we construct the reciprocal lattice of vectors $`\stackrel{}{K}`$. It is convenient to partition the space into the region $`S`$ inside the non-overlapping spheres of radius $`R`$ and the interstitial region $`I`$, Fig. 1.

We expand the charge density in Fourier series in region $`I`$ and in multipolar terms in region $`S`$ (dual representation ):

$`e\rho (\stackrel{}{R})|_{\stackrel{}{R}I}=e{\displaystyle \underset{\stackrel{}{K}}{}}\rho _I(\stackrel{}{K})e^{i\stackrel{}{K}\stackrel{}{R}},`$ (3)
$`e\rho (\stackrel{}{R})|_{\stackrel{}{R}S}=e\left({\displaystyle \underset{\mathrm{\Lambda }}{}}\rho _\mathrm{\Lambda }(r)S_\mathrm{\Lambda }(\widehat{r})Z{\displaystyle \frac{\delta (r)}{\sqrt{4\pi }r^2}}S_0\right).`$ (4)
Angular functions $`S_\mathrm{\Lambda }`$ are symmetry adapted functions  (SAFs) which are linear combinations of spherical harmonics (see Appendix A for details and definitions). The total density is an invariant of the space group under consideration. Therefore, in the expansion (4) we use only SAFs $`S_\mathrm{\Lambda }`$ of the full (or unit) $`A_{1g}`$ symmetry. Finally, $`\delta (r)`$ in (4) is the one dimensional (radial) delta-function, $`Z`$ is the atomic number and $`e<0`$ is the charge of electron. In the following we will use atomic units and the charge density is understood in units $`e=1`$. The density expansion (3,b) called also the dual representation  is characteristic of the LAPW method.  In fact, the dual representation is natural because in the interstitial region electron density is smooth whereas large charge oscillations near nuclei require a multipole analysis.

Using the linearity of Poisson’s equation we can extend smoothly the interstitial Fourier charge representation inside the spheres and subtract it out. We obtain

$`\rho (\stackrel{}{R})=\rho _I(\stackrel{}{R})+{\displaystyle \underset{S(\stackrel{}{n})}{}}\rho _{S(\stackrel{}{n})}^{}(\stackrel{}{R}),`$ (6)
$`\rho _{S(\stackrel{}{n})}^{}(\stackrel{}{R})=(\rho _{S(\stackrel{}{n})}(\stackrel{}{R})\rho _I(\stackrel{}{R}))\mathrm{\Theta }(\stackrel{}{R}S(\stackrel{}{n})),`$ (7)
where $`\rho _I(\stackrel{}{R})`$ is defined by (3) over the whole space, $`\rho _{S(\stackrel{}{n})}^{}`$ is the renormalized charge density inside a sphere $`S(\stackrel{}{n})`$ and $`\mathrm{\Theta }`$ is the unit step function. Using Eq. (1) and the plane wave expansion in SAFs centered at the site $`\stackrel{}{n}`$,

$`e^{i\stackrel{}{K}\stackrel{}{R}}=e^{i\stackrel{}{K}\stackrel{}{X}(\stackrel{}{n})}4\pi {\displaystyle \underset{\mathrm{\Lambda }}{}}i^lj_l(Kr(\stackrel{}{n}))S_\mathrm{\Lambda }(\widehat{K})S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n})),`$ (8)
we rewrite the multipolar density inside the sphere as

$`\rho _0^{}(r)=\rho _0(r)\sqrt{4\pi }\rho _I(\stackrel{}{K}=0)Z{\displaystyle \frac{\delta (r)}{\sqrt{4\pi }r^2}},`$ (10)
$`\rho _\mathrm{\Lambda }^{}(r)=\rho _\mathrm{\Lambda }(r)4\pi i^l{\displaystyle \underset{\stackrel{}{K}0}{}^{}}j_l(Kr)S_\mathrm{\Lambda }(\widehat{K})\rho _I(\stackrel{}{K}).`$ (11)
We shall refer to (6,b) with the multipole densities (10,b) as to the equivalent representation of the charge density. We will use it every time when we need to obtain the potential in the interstitial region because the Fourier series (3) is now extended to the whole crystal, that has some important advantages (see next subsection).

In principle, in Eq. (3) there is a term with $`\stackrel{}{K}=0`$ that corresponds to the homogeneous electron density distribution. However, it will not contribute to the resulting potential because it cancels with the positive homogeneous density distribution arising from nuclei (see also next subsection).

### B The solution in the interstitial region

The potential in the interstices has three components. The first contribution arises from the density $`\rho _I(\stackrel{}{R})`$ given by the Fourier series (3) throughout the whole crystal. The second and third terms are due to the monopole ($`l=0`$) and the high multipole ($`l1`$) moments of electron density inside the spheres. (We recall that here we are dealing with the equivalent charge density given by Eqs. (10,b).)

From the Fourier series of $`\rho _I(\stackrel{}{R})`$ we find

$`V_I(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{K}0}{}^{}}V^I(\stackrel{}{K})e^{i\stackrel{}{K}\stackrel{}{R}},`$ (12)
where

$`V^I(\stackrel{}{K})={\displaystyle \frac{4\pi }{K^2}}\rho _I(\stackrel{}{K}).`$ (13)
The easiness of the last expression owes to the fact that we extend the Fourier expansion of $`\rho _I(\stackrel{}{R})`$ to the whole crystal. A complication is that the density inside a sphere is renormalized according to Eqs. (10,b). Next we consider a sphere $`S(\stackrel{}{n})`$ and for a moment we will not include to the potential the contribution due to the charges outside the sphere. Thus we introduce a single sphere potential $`𝒰_S(\stackrel{}{r})=𝒰_{S(\stackrel{}{n})}(\stackrel{}{r}(\stackrel{}{n}))`$. Using the one-site expansion  in SAFs

$`{\displaystyle \frac{1}{|\stackrel{}{r}(\stackrel{}{n})\stackrel{}{r}^{}(\stackrel{}{n})|}}={\displaystyle \underset{\mathrm{\Lambda }}{}}v_{\mathrm{\Lambda }\mathrm{\Lambda }}(\stackrel{}{n}=\stackrel{}{n}^{},r,r^{})S_\mathrm{\Lambda }(\widehat{r}^{})S_\mathrm{\Lambda }(\widehat{r}),`$ (14)
where

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n}=\stackrel{}{n}^{},r,r^{})`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Omega }𝑑\mathrm{\Omega }^{}\frac{S_\mathrm{\Lambda }(\mathrm{\Omega })S_\mathrm{\Lambda }^{}(\mathrm{\Omega }^{})}{|\stackrel{}{r}(\stackrel{}{n})\stackrel{}{r}^{}(\stackrel{}{n})|}}`$ (15)
$`=`$ $`\delta _{\mathrm{\Lambda }\mathrm{\Lambda }^{}}{\displaystyle \frac{4\pi }{2l+1}}{\displaystyle \frac{r_<^l}{r_>^{l+1}}},`$ (16)
and $`r_<`$ ($`r_>`$) is the smaller (larger) of $`r`$ and $`r^{}`$, we obtain the solution inside the sphere ($`0<rR`$),

$`𝒰_S(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \frac{Q_0(r)}{r}}+Q_0^{}(r)`$ (17)
$`+`$ $`{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \frac{4\pi }{2l+1}}\left({\displaystyle \frac{Q_\mathrm{\Lambda }(r)}{r^{l+1}}}+r^lQ_\mathrm{\Lambda }^{}(r)\right)S_\mathrm{\Lambda }(\widehat{r}).`$ (18)
Here we have introduced the functions

$`Q_0(r)=q_0(r){\displaystyle \frac{4\pi r^3}{3}}\rho _I(\stackrel{}{K}=0),`$ (20)
$`Q_0^{}(r)=q_0^{}(r)2\pi (R^2r^2)\rho _I(\stackrel{}{K}=0),`$ (21)
with

$`q_0(r)=\sqrt{4\pi }{\displaystyle _0^r}\rho _0(r^{})r_{}^{}{}_{}{}^{2}𝑑r^{}Z,`$ (23)
$`q_0^{}(r)=\sqrt{4\pi }{\displaystyle _r^R}\rho _0(r^{})r^{}𝑑r^{},`$ (24)
for the spherically symmetric term $`\mathrm{\Lambda }=l=0`$, and

$`Q_\mathrm{\Lambda }(r)={\displaystyle _0^r}\rho _\mathrm{\Lambda }^{}(r^{})r_{}^{}{}_{}{}^{l+2}𝑑r^{},`$ (26)
$`Q_\mathrm{\Lambda }^{}(r)={\displaystyle _r^R}\rho _\mathrm{\Lambda }^{}(r^{})r_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}1}l}𝑑r^{},`$ (27)
for the other multipoles, $`\mathrm{\Lambda }0`$. Taking into account the explicit expression for $`\rho ^{}(r)`$, Eq. (11), we rewrite (26,b) as

$`Q_\mathrm{\Lambda }(r)`$ $`=`$ $`q_\mathrm{\Lambda }(r)4\pi i^lr^{l+2}{\displaystyle \underset{\stackrel{}{K}0}{}^{}}S_\mathrm{\Lambda }(\widehat{K})\rho _I(\stackrel{}{K}){\displaystyle \frac{j_{l+1}(Kr)}{K}},`$ (29)
$`Q_\mathrm{\Lambda }^{}(r)`$ $`=`$ $`q_\mathrm{\Lambda }^{}(r)+4\pi i^l{\displaystyle \underset{\stackrel{}{K}0}{}^{}}S_\mathrm{\Lambda }(\widehat{K})\rho _I(\stackrel{}{K})`$ (31)
$`\times `$ $`{\displaystyle \frac{1}{K}}\left({\displaystyle \frac{j_{l1}(KR)}{R^{l1}}}{\displaystyle \frac{j_{l1}(Kr)}{r^{l1}}}\right),`$ (32)
where

$`q_\mathrm{\Lambda }(r)`$ $`=`$ $`{\displaystyle _0^r}\rho _\mathrm{\Lambda }(r^{})r_{}^{}{}_{}{}^{l+2}𝑑r^{},`$ (33)
$`q_\mathrm{\Lambda }^{}(r)`$ $`=`$ $`{\displaystyle _r^R}\rho _\mathrm{\Lambda }(r^{})r_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}1}l}𝑑r^{}.`$ (34)
In order to obtain Eq. (LABEL:2.12a,b) we have performed integration of spherical Bessel functions using the properties 10.1.23, 10.1.24 of Ref. .

For the potential $`𝒰_S(\stackrel{}{r})`$ outside the sphere $`S(\stackrel{}{n})`$ ($`rR`$) we have

$`𝒰_S(\stackrel{}{r})={\displaystyle \frac{Q_0(R)}{r}}+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \frac{4\pi }{2l+1}}{\displaystyle \frac{Q_\mathrm{\Lambda }(R)}{r^{l+1}}}S_\mathrm{\Lambda }(\widehat{r}).`$ (35)
Now we are ready to compute the potential of the periodic arrangements of the spheres,

$`V_S(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{n}}{}}𝒰_S(\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})).`$ (36)
Since we seek for the solution in the interstitial region it is desirable to expand it in a Fourier series, in the spirit of the dual representation,

$`V_S(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{K}0}{}^{}}V^S(\stackrel{}{K})e^{i\stackrel{}{K}\stackrel{}{R}}.`$ (38)
We distinguish two contributions, $`V_0^S`$ and $`V_M^S`$, from the spherically symmetric charge distribution and from the other higher charge multipoles, respectively:

$`V^S(\stackrel{}{K})=V_0^S(\stackrel{}{K})+V_M^S(\stackrel{}{K}).`$ (39)
The monopole term ($`l=0`$) $`V_0^S`$ can not be computed directly and requires a special treatment. Here we employ the Ewald method  and obtain for the Fourier coefficients

$`V_0^S(\stackrel{}{K})={\displaystyle \frac{\pi }{v}}{\displaystyle \frac{4}{K^2}}e^{K^2/4G^2}Q_0(R),`$ (40)
where $`G`$ is the Ewald cut-off parameter and $`v`$ is the unit cell volume, $`v=V/N`$. The improved Ewald technique requires an additional summation in real space, see Ref. . Since we are concerned with the Fourier expansion in the interstitial region we can always take $`G`$ large enough and neglect the summation in real space. (At large $`G`$ the dimensionless parameter $`c=GR1`$ and the complementary error function decreases exponentially.) Then the Ewald expansion with the Fourier coefficients (40) converges absolutely in the interstitial region.

In the Ewald expansion, Eq. (40), and in the Fourier decomposition (12) we have omitted the terms with $`\stackrel{}{K}=0`$ which correspond to a homogeneous charge distribution of the electrons and the nuclei. Both these terms, considered separately, diverge. Their sum, however, as a consequence of electroneutrality gives the zero value of the charge density and therefore can be discarded.

The high multipole terms ($`l1`$) converge fast. Their Fourier coefficients are given by

$`V_M^S(\stackrel{}{K})={\displaystyle \frac{1}{v}}{\displaystyle _{cell}}e^{i\stackrel{}{K}\stackrel{}{R}}V_{S_M}(\stackrel{}{R})d^3R.`$ (41)
Here the integration is over any primitive unit cell of the crystal, and $`V_{S_M}(\stackrel{}{R})`$ is the potential (36) where only the high multipoles ($`l1`$) are retained in $`𝒰_S`$. Since the result of integration (41) does not depend on the choice of unit cell, it follows  that the coefficients $`V_M^S(\stackrel{}{K})`$ can be rewritten as

$`V_M^S(\stackrel{}{K})={\displaystyle \frac{1}{v}}{\displaystyle _{crystal}}e^{i\stackrel{}{K}\stackrel{}{R}}𝒰_{S_M}(\stackrel{}{R})d^3R.`$ (42)
Here the integration is taken over the whole crystal and $`𝒰_{S_M}(\stackrel{}{R})`$ is given by the last (high multipole) term on the right hand sides of (18) and (35) for $`0<|\stackrel{}{R}|R`$ and $`|\stackrel{}{R}|>R`$, respectively. Using Eq. 10.1.24 of Ref.  we find

$`V_M^S(\stackrel{}{K})`$ $`=`$ $`{\displaystyle \frac{(4\pi )^2}{v}}{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \frac{(i)^l}{2l+1}}S_\mathrm{\Lambda }(\widehat{K})`$ (45)
$`\times \left(A_\mathrm{\Lambda }(K)+Q_\mathrm{\Lambda }(R){\displaystyle \frac{j_{l1}(KR)}{KR^{l1}}}\right),`$
where

$`A_\mathrm{\Lambda }(K)={\displaystyle _0^R}\left({\displaystyle \frac{Q_\mathrm{\Lambda }(r)}{r^{l+1}}}+r^lQ_\mathrm{\Lambda }^{}(r)\right)j_l(Kr)r^2𝑑r.`$ (46)
(Note that at $`r0`$, $`Q_\mathrm{\Lambda }r^{l+2}`$, $`Q_\mathrm{\Lambda }^{}r^{1l}`$, $`j_l(Kr)r^l`$ and the integral (46) does not diverge.)

Collecting all contributions together, Eqs. (13), (40), (45), we obtain for the Coulomb potential in the interstitial region,

$`V(\stackrel{}{R}I)={\displaystyle \underset{\stackrel{}{K}0}{}^{}}V(\stackrel{}{K})e^{i\stackrel{}{K}\stackrel{}{R}},`$ (48)
where

$`V(\stackrel{}{K})=V^I(\stackrel{}{K})+V_0^S(\stackrel{}{K})+V_M^S(\stackrel{}{K}).`$ (49)
### C The solution inside the spheres

The potential in the form (48,b) can not be extended directly to the region $`S`$ inside the spheres. The main reason is that it requires an infinite number of plane waves to describe the Coulomb nucleus singularity. In order to obtain a well-behaved solution we expand the total potential $`V(\stackrel{}{R})`$ inside a sphere $`S(\stackrel{}{n})`$ ($`r(\stackrel{}{n})R`$) in a multipole series

$`V(\stackrel{}{r}(\stackrel{}{n}))=V_0(r)+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}V_\mathrm{\Lambda }(r)S_\mathrm{\Lambda }(\widehat{r}),`$ (50)
The potential $`V(\stackrel{}{r}(\stackrel{}{n}))`$ is defined uniquely by boundary conditions.  These are of the Dirichlet kind here, since by using the expansion (48) for $`V(\stackrel{}{R}I)`$ on the sphere surface, we arrive at

$`V_0(R)={\displaystyle \underset{\stackrel{}{K}0}{}^{}}j_0(KR)V(\stackrel{}{K}),`$ (52)
$`V_\mathrm{\Lambda }(R)=4\pi i^l{\displaystyle \underset{\stackrel{}{K}0}{}^{}}j_l(KR)S_\mathrm{\Lambda }(\widehat{K})V(\stackrel{}{K}).`$ (53)
Now inside the sphere $`S(\stackrel{}{n})`$ we have a Dirichlet problem.  We solve it by considering the spherical Green’s function expansion,  and rewrite the solution (Eqs. (3.126) and (3.114) of Ref. ) in the following form:

$`V(\stackrel{}{r}(\stackrel{}{n}))=U_S(\stackrel{}{r}(\stackrel{}{n}))+V_S^{out}(\stackrel{}{r}(\stackrel{}{n})),`$ (54)
The potential $`U_S(\stackrel{}{r}(\stackrel{}{n}))`$ is induced by the charge distribution inside the sphere $`S(\stackrel{}{n})`$. Clearly, it is independent of site $`\stackrel{}{n}`$ and we find

$`U_S(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \frac{q_0(r)}{r}}+q_0^{}(r)`$ (55)
$`+`$ $`{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \frac{4\pi }{2l+1}}\left({\displaystyle \frac{q_\mathrm{\Lambda }(r)}{r^{l+1}}}+r^lq_\mathrm{\Lambda }^{}(r)\right)S_\mathrm{\Lambda }(\widehat{r}).`$ (56)
Notice that the potential $`U_S`$ differs from $`𝒰_S`$, Eq. (35). In (56) we use $`q_\mathrm{\Lambda }`$, Eq. (33), and $`q_\mathrm{\Lambda }^{}`$, Eq. (34) (not $`Q_\mathrm{\Lambda }`$ and $`Q_\mathrm{\Lambda }^{}`$!). $`q_\mathrm{\Lambda }`$ and $`q_\mathrm{\Lambda }^{}`$ correspond to the initial charge density inside the sphere, Eq. (4), rather than the equivalent one, Eqs. (10,b). $`V_S^{out}(\stackrel{}{r})`$ is the potential due to the charge situated outside the sphere $`S(\stackrel{}{n})`$,

$`V_S^{out}(\stackrel{}{r})=V_0^{out}(R)+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}V_\mathrm{\Lambda }^{out}(R)\left({\displaystyle \frac{r}{R}}\right)^lS_\mathrm{\Lambda }(\widehat{r}).`$ (57)
Here the constants $`V_0^{out}V_0^{out}(R)`$ and $`V_\mathrm{\Lambda }^{out}V_\mathrm{\Lambda }^{out}(R)`$ are given by

$`V_0^{out}`$ $`=`$ $`V_0{\displaystyle \frac{q_0(R)}{R}},`$ (59)
$`V_\mathrm{\Lambda }^{out}`$ $`=`$ $`V_\mathrm{\Lambda }{\displaystyle \frac{4\pi }{2l+1}}{\displaystyle \frac{q_\mathrm{\Lambda }(R)}{R^{l+1}}}.`$ (60)
In the following we drop the argument $`R`$ in the boundary values $`V_0V_0(R)`$ and $`V_\mathrm{\Lambda }V_\mathrm{\Lambda }(R)`$. Notice that at $`\stackrel{}{r}=0`$ the potential due to the charges outside the sphere $`S(\stackrel{}{n})`$ is given only by $`V_0^{out}`$ since the other contributions ($`l0`$) reduce to zero.

Recalling that the Fourier components of the interstitial potential (48,b) consists of three parts, we distinguish the same three contributions on the sphere surface for $`V_0`$ and $`V_\mathrm{\Lambda }`$ in Eqs. (52,b):

$`V_\mathrm{\Lambda }=V_\mathrm{\Lambda }^I+V_\mathrm{\Lambda }^{S_0}+V_\mathrm{\Lambda }^{S_M}.`$ (61)
Here $`V_\mathrm{\Lambda }^I`$, $`V_\mathrm{\Lambda }^{S_0}`$ and $`V_\mathrm{\Lambda }^{S_M}`$ are obtained by using Eqs. (52,b) where instead of $`V(\stackrel{}{K})`$ we substitute $`V^I(\stackrel{}{K})`$, $`V_0^S(\stackrel{}{K})`$ and $`V_M^S(\stackrel{}{K})`$, respectively. Taking into account (59,b) we can separate the contributions to $`V_\mathrm{\Lambda }^{out}`$ into three parts,

$`V_\mathrm{\Lambda }^{out}=V_{\mathrm{\Lambda },I}^{out}+V_{\mathrm{\Lambda },S_0}^{out}+V_{\mathrm{\Lambda },S_M}^{out}.`$ (62)
Here $`V_{\mathrm{\Lambda },I}^{out}`$ is the contribution from the interstitial region, while $`V_{\mathrm{\Lambda },S_0}^{out}`$ and $`V_{\mathrm{\Lambda },S_M}^{out}`$ accounts for the contributions from the monopoles and the higher multipoles from the sites $`\stackrel{}{n}^{}\stackrel{}{n}`$, respectively. Explicitly, we obtain

$`V_{\mathrm{\Lambda },I}^{out}=V_\mathrm{\Lambda }^I,`$ (64)
$`V_{\mathrm{\Lambda },S_0}^{out}=\{\begin{array}{cc}V_0^{S_0}\frac{q_0(R)}{R},\hfill & \mathrm{\Lambda }=0,\hfill \\ V_\mathrm{\Lambda }^{S_0},\hfill & \mathrm{\Lambda }0,\hfill \end{array}`$ (67)
$`V_{\mathrm{\Lambda },S_M}^{out}=\{\begin{array}{cc}V_0^{S_M},\hfill & \mathrm{\Lambda }=0,\hfill \\ V_\mathrm{\Lambda }^{S_M}\frac{4\pi }{2l+1}\frac{q_\mathrm{\Lambda }(R)}{R^{l+1}},\hfill & \mathrm{\Lambda }0.\hfill \end{array}`$ (70)
The representation of $`V(\stackrel{}{R})`$ given by Eqs. (54) and (62) is convenient because one clearly distinguishes the intra-site part $`V_S`$ from the inter- site ones, $`V_{\mathrm{\Lambda },S_0}^{out}`$, $`V_{\mathrm{\Lambda },S_M}^{out}`$, and $`V_{\mathrm{\Lambda },I}^{out}`$.

### D Comparison with the two-center expansion of the Coulomb potential

It is instructive to compare the results of the previous subsection with the two-center expansion of the Coulomb potential,

$`{\displaystyle \frac{1}{|\stackrel{}{R}(\stackrel{}{n})\stackrel{}{R}^{}(\stackrel{}{n}^{})|}}`$ (71)
$`={\displaystyle \underset{\mathrm{\Lambda }\mathrm{\Lambda }^{}}{}}v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n}))S_\mathrm{\Lambda }^{}(\widehat{r}^{}(\stackrel{}{n}^{})),`$ (72)
where

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})={\displaystyle 𝑑\mathrm{\Omega }(\stackrel{}{n})𝑑\mathrm{\Omega }(\stackrel{}{n}^{})\frac{S_\mathrm{\Lambda }(\widehat{n})S_\mathrm{\Lambda }^{}(\widehat{n}^{})}{|\stackrel{}{R}(\stackrel{}{n})\stackrel{}{R}^{}(\stackrel{}{n}^{})|}}.`$ (73)
(74)
The potential on the sphere surface at a site $`\stackrel{}{n}`$ is found by integrating (72) over $`\stackrel{}{r}^{}(\stackrel{}{n}^{})`$ and summing over the spheres $`\stackrel{}{n}^{}(\stackrel{}{n})`$:

$`V_S^{out}(\mathrm{\Omega })=V_0^{out}+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}(V_{\mathrm{\Lambda },S_0}^{out}+V_{\mathrm{\Lambda },S_M}^{out})S_\mathrm{\Lambda }(\mathrm{\Omega }),`$ (75)
where (compare with Eqs. (67,c))

$`V_{\mathrm{\Lambda },S_0}^{out}={\displaystyle \underset{\stackrel{}{n}^{}\stackrel{}{n}}{}^{}}v_{\mathrm{\Lambda }\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}^{},R),`$ (77)
$`V_{\mathrm{\Lambda },S_M}^{out}={\displaystyle \underset{\stackrel{}{n}^{}\stackrel{}{n}}{}^{}}{\displaystyle \underset{\mathrm{\Lambda }^{}0}{}^{}}v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{},R),`$ (78)
and

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{},R)={\displaystyle _0^R}v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};R,r^{})\rho _\mathrm{\Lambda }^{}^{}(r^{})r^2𝑑r^{}.`$ (79)
Since

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{}){\displaystyle \frac{(r)^l(r^{})^l^{}}{|\stackrel{}{X}(\stackrel{}{n})\stackrel{}{X}(\stackrel{}{n}^{})|^{l+l^{}+1}}},`$ (80)
the sums on $`\stackrel{}{n}^{}`$, Eqs. (77,b), converge fast (except for the case $`l=l^{}=0`$ which is irrelevant here). Practically, it is enough to consider summations over a few neighboring shells.

In many cases (for the cubic symmetry, for example) the whole sum (78) is small in comparison with (77). Eq. (77) then corresponds to the potential calculated with spherical symmetric charge distributions for the neighbors ($`i.e.`$ $`\mathrm{\Lambda }^{}=l^{}=0`$) and a homogeneous charge distribution in the interstices, Eqs. (10), (20). It follows from (80) that $`v_{\mathrm{\Lambda }0}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})`$ does not depend on $`r^{}`$ which is a dummy argument and

$`v_{\mathrm{\Lambda }\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}^{},R)=B_\mathrm{\Lambda }(\stackrel{}{n}^{})Q_0(R){\displaystyle \frac{1}{\sqrt{4\pi }}}.`$ (81)
Here $`Q_0(R)`$ is the effective charge inside the sphere given by Eq. (20), and $`B_\mathrm{\Lambda }(\stackrel{}{n}^{})=v_{\mathrm{\Lambda }0}(\stackrel{}{n},\stackrel{}{n}^{};R,r^{})`$. Combining (77) and (81) with (57), we introduce the potential

$`V^{CEF}(\stackrel{}{r})={\displaystyle \frac{Q_0(R)}{\sqrt{4\pi }}}{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}B_\mathrm{\Lambda }\left({\displaystyle \frac{r}{R}}\right)^lS_\mathrm{\Lambda }(\widehat{r}),`$ (83)
where

$`B_\mathrm{\Lambda }={\displaystyle \underset{\stackrel{}{n}^{}}{}^{}}B_\mathrm{\Lambda }(\stackrel{}{n}^{}).`$ (84)
The potential $`V^{CEF}`$ represents the crystal electric field potential (CEF) for localized electrons at site $`\stackrel{}{n}`$ if there are no multipolar couplings with conduction electrons and the electron density in the interstitial region is homogeneous, $`V^I(\stackrel{}{K})=0`$. Examples of calculations of $`B_\mathrm{\Lambda }`$ for a face centered cubic lattice (fcc) will be given in the next section. In the nearest neighbor approximation one has $`B_\mathrm{\Lambda }=N_{nn}B_\mathrm{\Lambda }(\stackrel{}{n}^{})`$, where $`\stackrel{}{n}^{}`$ refers to any of the nearest neighbors and $`N_{nn}`$ is the number of nearest neighbors. Such CEF potential has been considered in Refs. . A more sophisticated CEF for localized and itinerant electrons will be discussed in Sec. IV.

## III Application to the face centered cubic lattice

In this section we apply the technique of Sec. II to test the method and to calculate some important structural constants. As a Bravais lattice we take the face centered cubic (fcc) lattice.

### A  Periodic Coulomb potential of the monopole ($`l=0`$) charge density

We first consider the potential of the unit ($`q=1`$) point charges situated on a face centered cubic lattice. For simplicity we take the cubic lattice constant $`a=1`$. (The volume of the unit cell is $`v=1/4`$.) We consider touching spheres on the fcc lattice with the close contact radius $`R=\sqrt{2}/4`$. As we know from Sec. II, the crystal potential is defined only if the total charge of the unit cell is zero. Therefore, complementary to the positive point charges we must introduce a compensating negative charge density distributed homogeneously in the interstitial region. The total charge confined in the interstitial region of the unit cell must be $`q_I=1`$ with the density $`\rho _I=1/v_I=4/(1\pi /3\sqrt{2})`$.

Using the method of Sec. II we first expand the potential of the interstitial region in Fourier series and then extend it smoothly to the whole crystal, Eqs. (2.17a,b). There are no high multipole moments inside a sphere and from Eq. (45) we find $`V_M^S(\stackrel{}{K})=0`$. Due to electroneutrality the terms with $`\stackrel{}{K}=0`$ are omitted. Since in the interstitial region there is only one component $`\rho _I(\stackrel{}{K}=0)`$ which is different from zero, this results in $`V^I(\stackrel{}{K}0)=0`$. Therefore, the potential in the interstices is given by (48) where the only nonzero contribution to the Fourier component $`V(\stackrel{}{K})`$ is the Ewald contribution $`V_0^S(\stackrel{}{K})`$, Eq. (40). Notice that for the Ewald expansion we use an effective point charge $`Q_0Q_{eff}=q\rho _I4\pi R^3/3=1+\pi /(3\sqrt{2}\pi )`$. The renormalization of charges is the price that we have to pay for the extension of the Fourier expansion inside the spheres, Sec. II.

The potential inside a sphere $`S(\stackrel{}{n})`$ then is given by

$`V(\stackrel{}{r})={\displaystyle \frac{1}{r}}+V_0^{out}+V_{04}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r})`$ $`+`$ $`V_{06}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r})`$ (86)
$`+\mathrm{}.`$
Here we introduce the cubic harmonics $`K_4(\mathrm{\Omega })`$ and $`K_6(\mathrm{\Omega })`$ which are SAFs $`S_{l=4,A1g}(\mathrm{\Omega })`$ and $`S_{l=6,A1g}(\mathrm{\Omega })`$, correspondingly,  and

$`V_0^{out}=V_0{\displaystyle \frac{1}{R}}.`$ (87)
It is interesting to remark that although we have started with the model of point charges on the fcc lattice, the total crystal potential (86) acquires contributions of the high multipole orders, $`l=4,6,8,\mathrm{}`$. The constants $`V_0`$, $`V_{04}`$, $`V_{06}`$ of the potential are found from the boundary conditions on the sphere surface, Eqs. (52,b). The potential energy of the charges (fcc lattice) is given by

$`E_{fcc}/N=V_0^{out}+{\displaystyle \frac{1}{2}}V_I=3.9885\text{ a.u.},`$ (88)
where $`V_I`$ is the average potential of the interstitial region $`I`$,

$`V_I={\displaystyle \frac{1}{v_I}}{\displaystyle _I}V(\stackrel{}{R})𝑑\stackrel{}{R}={\displaystyle \frac{v}{v_I}}{\displaystyle \underset{\stackrel{}{K}0}{}^{}}V(\stackrel{}{K})𝒪(\stackrel{}{K}).`$ (89)
The overlap integral $`𝒪(\stackrel{}{K})`$ reads

$`𝒪(\stackrel{}{K})=\delta _\stackrel{}{K}{\displaystyle \frac{4\pi R^2}{v}}{\displaystyle \frac{j_1(KR)}{K}}.`$ (90)
Thus $`E_{fcc}/N`$, Eq. (88), has the meaning of the Madelung constant for the fcc lattice.

The values $`V_0`$, $`V_{04}`$, $`V_{06}`$ and $`V_I`$ are proportional to the effective charge $`Q_{eff}1`$. It is convenient to consider the normalized to the unit charge quantities $`\stackrel{~}{V}_0=V_0/Q_{eff}`$, $`\stackrel{~}{V}_{04}=V_{04}/Q_{eff}`$, $`\stackrel{~}{V}_{06}=V_{06}/Q_{eff}`$ and $`\stackrel{~}{V}_I=V_I/Q_{eff}`$. We study the convergence of $`\stackrel{~}{V}_0`$, $`\stackrel{~}{V}_{04}`$, $`\stackrel{~}{V}_{06}`$ and $`\stackrel{~}{V}_I`$ by employing the Ewald method and the technique of Weinert, Ref. . The results are quoted in Tables I and II where the same number of reciprocal vectors $`|\stackrel{}{K}|<K_{max}`$ is used for both approaches.

(Though we did not optimize the convergence of Fourier series by choosing the best $`G`$ and $`n`$ for the Ewald and Weinert approaches.) We observe that $`\stackrel{~}{V}_{06}`$ and $`\stackrel{~}{V}_{04}`$ converge faster than $`\stackrel{~}{V}_0`$ and $`\stackrel{~}{V}_I`$.

The Ewald method estimates better the monopole terms $`\stackrel{~}{V}_0`$ and $`\stackrel{~}{V}_I`$, Table I. In this case it is simple and stable, while Weinert’ procedure requires evaluation of spherical Bessel functions of higher order for each $`\stackrel{}{K}`$.

Both methods give almost the same results for $`\stackrel{~}{V}_{04}`$ and $`\stackrel{~}{V}_{06}`$, Table II. $`\stackrel{~}{V}_{04}`$ and $`\stackrel{~}{V}_{06}`$ can be calculated independently by means of the two-center expansion of Coulomb interaction, which we have considered in Sec. II D. In that case we first calculate $`v_{04}(\stackrel{}{n},\stackrel{}{n}^{})`$ and $`v_{06}(\stackrel{}{n},\stackrel{}{n}^{})`$, Eq. (74), where $`r=r^{}=R`$, with $`S_0=1/\sqrt{4\pi }`$ on a first site $`\stackrel{}{n}`$, and the cubic harmonic (SAF) $`K_4`$ or $`K_6`$ on a second site $`\stackrel{}{n}^{}`$. Since the integrals $`v_{04}(\stackrel{}{n},\stackrel{}{n}^{})`$ and $`v_{06}(\stackrel{}{n},\stackrel{}{n}^{})`$ decrease very fast with the distance $`d`$ between the two sites ($`1/d^5`$ for $`K_4`$ and $`1/d^7`$ for $`K_6`$), we obtain $`\stackrel{~}{V}_{04}`$ and $`\stackrel{~}{V}_{06}`$ by summing $`v_{04}(\stackrel{}{n},\stackrel{}{n}^{})`$ and $`v_{06}(\stackrel{}{n},\stackrel{}{n}^{})`$ over all neighbors $`\stackrel{}{n}^{}`$ on the fcc lattice, i.e.

$`\stackrel{~}{V}_{0t}={\displaystyle \frac{1}{\sqrt{4\pi }}}{\displaystyle \underset{\stackrel{}{n}^{}}{}^{}}v_{0t}(\stackrel{}{n},\stackrel{}{n}^{}),`$ (91)
where $`t=4`$ or 6. The results of such calculations are shown in Tables III and IV.

The convergent values $`\stackrel{~}{V}_{04}`$ and $`\stackrel{~}{V}_{06}`$ perfectly match those obtained by the techniques of Ewald and Weinert, Table II. However, the Ewald and Weinert’ approaches require less computer time and are more efficient than the two-center expansion. Finally, we would like to remark that the calculated in Table II exact values $`\stackrel{~}{V}_{04}`$ and $`\stackrel{~}{V}_{06}`$ can be used for cubic crystal electric field, Eq. (83,b). For an fcc crystal with the cubic lattice constant $`a1`$ we have $`B_4=\sqrt{4\pi }\stackrel{~}{V}_{04}/a`$ and $`B_6=\sqrt{4\pi }\stackrel{~}{V}_{06}/a`$.

Although we have analyzed both methods at the close contact radius $`R=\sqrt{2}/4`$ it is also possible to assess the validity of the two methods by decreasing the sphere radius $`rR`$. Both methods are designed to describe the potential in the interstitial region. They can not approximate the potential correctly in the whole space. Therefore, both of them are expected to fail at a certain small radius, but their departure from the exact solution is an indication of their accuracy. The results for $`\stackrel{~}{V}_{04}(r)`$ and $`\stackrel{~}{V}_{06}(r)`$ ($`K_{max}R`$=90.73, $`G`$=40, $`n`$=84) at $`r<R`$ are shown in Fig. 2 and 3.

Weinert’ procedure gives more smooth values for $`r`$ close to $`R`$, but the Ewald expansion is more stable in the whole range. The calculation of the monopole term $`V_0^{out}`$ (or $`V_0`$) is of special interest. By decreasing the sphere radius $`r`$ we change the density of the compensating negative charge distribution. As a result, both $`V_0(r)`$ and $`V_0^{out}(r)`$, Eq. (87), are functions of $`r`$. From the last value however it is possible to reconstruct $`V_0^{out}(R)`$ at the close contact radius $`R`$. On can easily obtain the result from electrostatic considerations:

$`V_0^{out}(R)|_r={\displaystyle \frac{V_0^{out}(r)+2\pi (R^2r^2)\rho }{1\frac{4}{3}\pi (R^3r^3)\rho }},`$ (93)
where

$`\rho =\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{4\pi }{3}}r^3\right)^1,`$ (94)
In principle, $`V_0^{out}(R)`$ should be independent of $`r`$, but in practice it exhibits such dependence as indicated on the left hand side of (93). This is because the estimation (93) depends on the accuracy of the calculated value $`V_0^{out}(r)`$ and as such, it is a function of $`r`$. We find it convenient to study $`V_0^{out}(R)|_r`$ which is expected to be constant rather than the initial value $`V_0^{out}(r)`$. In order to compare the two techniques we plot the function $`\mathrm{}(r)=(V_0^{out}|_rV_0^{out}|_{ex})/V_0^{out}|_{ex}`$ in Fig. 4. ($`V_0^{out}|_{ex}=5.5612`$ is the best value for $`V_0^{out}`$ calculated from data of Table I, last row.) As one clearly sees from Fig. 4, the Ewald expansion performs better than the approach of Weinert.

### B Periodic Coulomb potentials of the $`l=4`$ and $`l=6`$ cubic charge density

Here we illustrate how the proposed method works for high order multipoles with $`l=4`$ and $`l=6`$. We consider two cases. In the first case the charge distribution on the touching sphere surfaces is given by

$`\rho _4(\stackrel{}{r})={\displaystyle \frac{\delta (rR)}{r^2}}K_4(\widehat{r}),`$ (96)
In the second case it is

$`\rho _6(\stackrel{}{r})={\displaystyle \frac{\delta (rR)}{r^2}}K_6(\widehat{r}).`$ (97)
Here $`\delta (r^{})`$ is the one-dimensional delta-function and $`K_4`$, $`K_6`$ are cubic harmonics (SAFs).  Now it is not necessary to introduce a compensating charge in the interstitial region since the total charge of the unit cell given by integration of (96) or (97) over the polar angles yields zero. According to the treatment of Sec. II we start with the potential of a single sphere, say $`S(\stackrel{}{n}=0)`$. The potentials then are

$`U_S^4(\stackrel{}{r})={\displaystyle \frac{4\pi }{9}}{\displaystyle \frac{1}{R}}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r}),`$ (99)
$`U_S^6(\stackrel{}{r})={\displaystyle \frac{4\pi }{13}}{\displaystyle \frac{1}{R}}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r}),`$ (100)
for $`rR`$, and

$`U_S^4(\stackrel{}{r})={\displaystyle \frac{4\pi }{9}}{\displaystyle \frac{1}{R}}\left({\displaystyle \frac{R}{r}}\right)^5K_4(\widehat{r}),`$ (102)
$`U_S^6(\stackrel{}{r})={\displaystyle \frac{4\pi }{13}}{\displaystyle \frac{1}{R}}\left({\displaystyle \frac{R}{r}}\right)^7K_6(\widehat{r}),`$ (103)
for $`rR`$. Notice that the potentials $`U_S^4`$ and $`U_S^6`$ have discontinuous derivatives at $`R`$. (Here we do not make any distinction between $`U(r)`$ and $`𝒰(r)`$, Sec. II, since there is no renormalization of the charge density inside a sphere.) Using the method of Sec. II we find from Eqs. (45,b) the Fourier transforms of the potentials of the periodic arrangements of such charged surfaces on the fcc lattice,

$`V_4(\stackrel{}{K})={\displaystyle \frac{(4\pi )^2}{9v}}{\displaystyle \frac{R}{K}}(j_3(KR)+j_5(KR))K_4(\stackrel{}{K}),`$ (105)
$`V_6(\stackrel{}{K})={\displaystyle \frac{(4\pi )^2}{13v}}{\displaystyle \frac{R}{K}}(j_5(KR)+j_7(KR))K_6(\stackrel{}{K}).`$ (106)
Notice, that we arrive at the same results if first we Fourier transform the densities,

$`\rho _4(\stackrel{}{K})={\displaystyle \frac{4\pi }{v}}j_4(KR)K_4(\stackrel{}{K}),`$ (108)
$`\rho _6(\stackrel{}{K})={\displaystyle \frac{4\pi }{v}}j_6(KR)K_6(\stackrel{}{K}),`$ (109)
and then calculate the potentials by means of

$`V_4(\stackrel{}{K})={\displaystyle \frac{4\pi }{K^2}}\rho _4(\stackrel{}{K}),`$ (111)
$`V_6(\stackrel{}{K})={\displaystyle \frac{4\pi }{K^2}}\rho _6(\stackrel{}{K}).`$ (112)
(In oder to establish the equality with (105,b) we used the property 10.1.21 of Bessel functions, Ref. .) As we have many times mentioned before (Sec. II) the Fourier expansion is valid only in the interstitial region.

Inside a sphere $`S(\stackrel{}{n})`$ the first potential is given by

$`V_4(\stackrel{}{r})=V_{40}+V_{44}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r})+V_{46}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r})+\mathrm{}.`$ (113)
(114)
For the second potential, accordingly, we have

$`V_6(\stackrel{}{r})=V_{60}+V_{64}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r})+V_{66}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r})+\mathrm{}.`$ (115)
(116)
The constants $`V_{40}`$, $`V_{60}`$, and $`V_{44}`$, $`V_{46}`$, $`V_{64}`$, $`V_{66}`$ are found from the boundary conditions, Eqs. (52) and (53), where $`V(\stackrel{}{K})`$ is given by $`V_4(\stackrel{}{K})`$, Eq. (105), for the first potential and by $`V_6(\stackrel{}{K})`$, Eq. (106), for the second potential. (Compare these results with the potentials (99) and (100) for a single sphere.)

In order to distinguish the potential due to the multipole moment of the sphere $`S(\stackrel{}{n})`$ from the potential of the rest of the crystal, we subtract from $`V_4`$ and $`V_6`$ the contribution induced by the multipole moment of the sphere (see Eqs. (64-c)) and introduce the parameters

$`V_{44}^{out}=V_{44}{\displaystyle \frac{4\pi }{9}}{\displaystyle \frac{1}{R}},`$ (118)
$`V_{66}^{out}=V_{66}{\displaystyle \frac{4\pi }{13}}{\displaystyle \frac{1}{R}}.`$ (119)
We have calculated $`V_{44}^{out}`$, $`V_{66}^{out}`$, $`V_{46}`$ and $`V_{60}`$ and then compared them with those obtained by the method of Ref. , see Tables V and VI.

Alternatively, we can compute $`V_{44}^{out}`$ and $`V_{66}^{out}`$ by performing the two-center multipole expansions (72), (74). To do it, we first calculate integrals $`v_{44}(\stackrel{}{n},\stackrel{}{n}^{})`$ and $`v_{66}(\stackrel{}{n},\stackrel{}{n}^{})`$, Eq. (74), between the functions $`K_4(\stackrel{}{n})`$ and $`K_4(\stackrel{}{n}^{})`$, or between $`K_6(\stackrel{}{n})`$ and $`K_6(\stackrel{}{n}^{})`$ at sites $`\stackrel{}{n}`$ and $`\stackrel{}{n}^{}`$. By summing over $`\stackrel{}{n}^{}`$ on the fcc lattice we estimate $`V_{44}^{out}`$ and $`V_{66}^{out}`$, i.e.

$`V_{tt}^{out}={\displaystyle \underset{\stackrel{}{n}^{}}{}^{}}v_{tt}(\stackrel{}{n},\stackrel{}{n}^{}),`$ (120)
where $`t=4`$ or 6. Clearly, the accuracy depends on the number of sites $`\stackrel{}{n}^{}`$ included in the summations. The results of these calculations are given in Table VII. (The representative coordinates of five nearest shells as well as the number of atoms of the shells are quoted in Tab. III.)

From the Tables V-VII it follows that the method of Ref.  gives a better convergence despite the fact that the expressions (108,b), (111,b) for the Fourier coefficients are exact. Interestingly, the exact expressions formally correspond to $`n_4=n_6=1`$ in Eq. 28 of Ref. , while the better convergence is achieved for large positive $`n_4`$ and $`n_6`$ implied by the technique of Ref. . We ascribe this property to the fact that the potentials $`V_4`$ and $`V_6`$ have discontinues derivatives at the sphere boundaries. Therefore, the exact expansions (108,b), (111,b) try to reproduce these cusps while the pseudo-charge density of Ref.  has $`n1`$ continuous derivatives thereby avoiding the problem. Thus, our test calculations for the densities (96) and (97) indicate that Weinert’ procedure gives a better convergence.

We can also generalize the present consideration for arbitrary multipole sphere moments $`Q_4(R)`$ and $`Q_6(R)`$ given by Eq. (26) on an fcc lattice when its cubic lattice constant $`a1`$. For the former case we obtain that the potential inside a sphere due to the other spheres’ multipoles $`Q_4`$ of the fcc crystal is

$`V_4^{out}(\stackrel{}{r})=B_{40}+B_{44}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r})+B_{46}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r})+\mathrm{},`$ (121)
(122)
where $`B_{40}=k_4V_{40}`$, $`B_{44}=k_4V_{44}^{out}`$, $`B_{46}=k_4V_{46}`$, and the coefficient of proportionality for touching spheres ($`R=\sqrt{2}a/4`$) is given by

$`k_4={\displaystyle \frac{64}{a^5}}Q_4(R).`$ (123)
Analogously, the expression for the potential $`V_6^{out}`$ due to $`Q_6(R)`$ is obtained as

$`V_6^{out}(\stackrel{}{r})=B_{60}+B_{64}\left({\displaystyle \frac{r}{R}}\right)^4K_4(\widehat{r})+B_{66}\left({\displaystyle \frac{r}{R}}\right)^6K_6(\widehat{r})+\mathrm{},`$ (124)
(125)
where $`B_{60}=k_6V_{40}`$, $`B_{64}=k_6V_{64}`$, $`B_{66}=k_6V_{66}^{out}`$, and

$`k_6={\displaystyle \frac{512}{a^7}}Q_6(R).`$ (126)
Notice that here we can immediately determine the potentials using exact values $`V_{40}=\stackrel{~}{V}_{04}`$, $`V_{44}`$, $`V_{46}`$, $`V_{60}=\stackrel{~}{V}_{06}`$, $`V_{64}`$ and $`V_{66}`$ quoted in Tables V and VI. These parameters as well as those calculated in Sec. III.A are in fact important structural constants of the fcc lattice.

## IV Hartree-Fock-Roothaan method: Calculations of the direct matrix elements

The Hartree-Fock operator for the electronic system is defined as

$`={\displaystyle \underset{d}{}}\left(_d+𝒥_d𝒦_d\right),`$ (127)
where $`_d`$ is the Hamiltonian operator for a $`d`$th electron moving in the field of nuclei alone,

$`_d={\displaystyle \frac{\stackrel{}{}_d^2}{2}}{\displaystyle \underset{\stackrel{}{n}}{}}{\displaystyle \frac{Z}{|\stackrel{}{R}_d\stackrel{}{X}(\stackrel{}{n})|}},`$ (128)
$`𝒥_d`$ and $`𝒦_d`$ are the direct Coulomb and exchange operators, respectively. In order to obtain the best trial function of the $`d`$th electron, $`\psi _d(x)=x|d`$, one considers basis functions $`\chi _a(x)=x|a`$, and expands $`\psi _d`$ in terms of $`\chi _a`$,

$`\psi _d(x)={\displaystyle \gamma _{da}\chi _a(x)},`$ (129)
where $`\gamma _{da}`$ are the coefficients of the expansion, and the index $`x`$ stands for the coordinates $`\stackrel{}{R}`$ and the spin projection $`s_z=\pm 1/2`$ of the electron, $`x(\stackrel{}{R},s_z)`$. For the matrix elements of $`𝒥_d`$ and $`𝒦_d`$ one has

$`a|𝒥_d|b={\displaystyle \underset{c}{}^{}}{\displaystyle \chi _a^{}(x_1)\chi _b(x_1)\frac{1}{R_{12}}\psi _c^{}(x_2)\psi _c(x_2)𝑑x_1𝑑x_2},`$ (131)
(132)
$`a|𝒦_d|b={\displaystyle \underset{c}{}^{}}{\displaystyle \chi _a^{}(x_1)\psi _c(x_1)\frac{1}{R_{12}}\chi _b(x_2)\psi _c^{}(x_2)𝑑x_1𝑑x_2}.`$ (133)
(134)
In Eqs. (132,b) the summation is understood over all occupied electron states except $`d`$. In the following we consider the restricted HF method. Usually in the HF method  the wave function $`\psi _d(x)=\varphi (\stackrel{}{R})\zeta (s_z)`$ is given by the coordinate (orbital) part $`\varphi (\stackrel{}{R})`$ and the spinor part $`\zeta (s_z)`$. However, it is well known that the spin-orbit coupling mixes the two components $`\zeta _+=\zeta (1/2)`$ and $`\zeta _{}=\zeta (1/2)`$. (The spin-orbit coupling is especially important for core shells.) Therefore, here we consider more general spin-orbitals,

$`\psi _d(x)=\varphi _+(\stackrel{}{R})\zeta _++\varphi _{}(\stackrel{}{R})\zeta _{}.`$ (135)
Using the time reversal symmetry $`𝒯`$,  one can construct from $`\psi _d`$ the time reversed state $`\psi _d^{}=𝒯\psi _d`$ (see Appendix C). According to the Kramers theorem $`\psi _d`$ and $`\psi _d^{}`$ have the same energy $`ϵ_d`$ and are orthogonal. Therefore, such treatment is analogous to the conventional restricted HF method  with the double occupancy of $`ϵ_d`$. However, in the following we will treat these two states separately, as different components of a doubled valued irreducible representation. Therefore, although we work within the restricted HF method, there are no factor 2 in front of $`_d`$ and $`𝒥_d`$ in Eq. (127). Notice, that as a consequence of (135), the integrals (132,b) include summations over two spin components.

We regroup $``$ as

$`=^{Coul}^{exc},`$ (136)
where

$`^{Coul}={\displaystyle \underset{d}{}}(_d+𝒥_d),`$ (137)
and

$`^{exc}={\displaystyle \underset{d}{}}𝒦_d.`$ (138)
The matrix elements of $`^{Coul}`$ can be written in the following form:

$`a|^{Coul}|b={\displaystyle \chi _a^{}(x)\left(\frac{1}{2}\stackrel{}{}^2+V^d(\stackrel{}{R})\right)\chi _b(x)𝑑x},`$ (139)
where the Hartree (electrostatic) potential $`V^d`$ acting on the $`d`$th electron reads

$`V^d(\stackrel{}{R})={\displaystyle \underset{c}{}^{}}{\displaystyle \frac{\psi _c^{}(x^{})\psi _c(x^{})}{|\stackrel{}{R}\stackrel{}{R}^{}|}𝑑x^{}}{\displaystyle \underset{\stackrel{}{n}}{}}{\displaystyle \frac{Z}{|\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})|}}.`$ (140)
(141)
It is convenient to rewrite the potential as

$`V^d(\stackrel{}{R})=V(\stackrel{}{R})𝒱_{}^{}{}_{}{}^{d}(\stackrel{}{R}),`$ (142)
where $`V(\stackrel{}{R})`$ is the potential of all electrons and nuclei considered in section II,

$`V(\stackrel{}{R})={\displaystyle \frac{\rho (\stackrel{}{R}^{})}{|\stackrel{}{R}\stackrel{}{R}^{}|}𝑑\stackrel{}{R}^{}},`$ (143)
while $`𝒱_{}^{}{}_{}{}^{d}`$ is the potential created by the $`d`$th electron alone,

$`𝒱_{}^{}{}_{}{}^{d}(\stackrel{}{R})={\displaystyle \frac{\psi _d^{}(x^{})\psi _d(x^{})}{|\stackrel{}{R}\stackrel{}{R}^{}|}𝑑x^{}}.`$ (144)
The charge density $`\rho (\stackrel{}{R})`$ in (143),

$`\rho (\stackrel{}{R})=\rho _{el}(\stackrel{}{R})+\rho _n(\stackrel{}{R}),`$ (146)
comprises the point charges of nuclei,

$`\rho _n(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{n}}{}}Z\delta (\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})),`$ (147)
and the total electronic charge density,

$`\rho _{el}(\stackrel{}{R})={\displaystyle \underset{c}{}}\psi _c^{}(x)\psi _c(x).`$ (148)
(Again, in (148) the summation is implied over two spin components.) In solids we are dealing with the two types of electron states, which are extended states of valence electrons (with the total charge density $`\rho _{val}(\stackrel{}{R})`$) and localized states of core electrons ($`\rho _{core}(\stackrel{}{r}(\stackrel{}{n}))`$), and

$`\rho _{el}(\stackrel{}{R})=\rho _{val}(\stackrel{}{R})+{\displaystyle \underset{\stackrel{}{n}}{}}\rho _{core}(\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})).`$ (149)
The Hartree-Fock-Roothaan method implies a self-consistent field procedure.  In particular, for itinerant states ($`d=(\stackrel{}{k},\alpha )`$) one solves the secular equation

$`{\displaystyle \underset{a}{}}_{ba}(\stackrel{}{k},\alpha )\gamma _a(\stackrel{}{k},\alpha )=E(\stackrel{}{k},\alpha ){\displaystyle \underset{a}{}}𝒪_{ba}(\stackrel{}{k})\gamma _a(\stackrel{}{k},\alpha ).`$ (150)
(151)
Here $`\stackrel{}{k}`$ is the wave vector, $`\alpha `$ is the band index, $`𝒪`$ is the overlap matrix and $`E(\stackrel{}{k},\alpha )`$ is the corresponding HF energy.  The HF operator $``$ depends implicitly on the solution (via the coefficients $`\gamma _a`$). Below in this section we calculate the matrix elements of $`^{Coul}`$, while in Sec. V we consider the matrix elements of $`^{exc}`$.

### A Multipole expansion of electron density

In order to calculate the matrix elements of $`^{Coul}`$ for a $`d`$th electron from Eqs. (139) we have to know the electrostatic potential $`V^d`$, Eq. (141) or (142). We have already thoroughly studied this problem in Sec. II and know how to proceed starting with the dual representation (3,b) of charge density. The only problem is to construct the charge density of the itinerant and the localized electrons. In this subsection we do it explicitly for the core (A1) and the valence (A2) electrons.

A1. The density of the inner closed shells at site $`\stackrel{}{n}`$ reads

$`\rho _{core}(\stackrel{}{r}(\stackrel{}{n}))={\displaystyle \underset{\tau }{}}n_\tau \rho _\tau (\stackrel{}{r}(\stackrel{}{n})),`$ (153)
where

$`\rho _\tau (\stackrel{}{r})=|\stackrel{}{r},s_z|\tau |^2=|x|\tau |^2.`$ (154)
Here $`x|\tau `$ stands for the wave functions of localized electrons, and $`n_\tau =1`$ if $`ϵ(\tau )<E_F`$ and zero otherwise ($`E_F`$ is the Fermi energy). In (153) summation over $`s_z=\pm 1/2`$ is implied. Here we will not consider the case of partially filled core shells. (Such situation occurs in lanthanides with localized 4$`f`$ electrons or in transition elements with $`d`$ electrons.) Furthermore, we assume that all core electrons at site $`\stackrel{}{n}`$ are confined inside the sphere $`S(\stackrel{}{n})`$ so that the core wave functions with their first derivatives fall to zero at $`R`$.

The inner closed shells are classified according to the principal quantum number $`n`$, the total angular momentum (including spin) $`j`$ and the orbital number $`l`$. In the case of the spherical symmetry these $`(2j+1)`$ electronic states belong to a doubled valued irreducible representation $`D_j`$ of the 3-dimensional rotation group $`SO(3)`$.  In the presence of crystal environment these levels (except $`s`$) in general are split into doubled valued irreducible representations $`(\mathrm{\Gamma },\nu )`$ of a crystal double group,

$`D_j{\displaystyle \mathrm{\Gamma }(\nu )},`$ (155)
where $`\nu `$ is used to label representations which occur more than once. We shall classify these core electronic states according to the atomic indices $`n`$, $`l`$, $`j`$, and $`\mathrm{\Gamma }`$, $`\nu `$, $`k`$ ($`k`$ labels the rows of $`\mathrm{\Gamma }`$) i.e. $`\tau (n,l,j;\mathrm{\Gamma },\nu ,k)`$, and

$`x|\tau =_\tau (r)\widehat{r},s_z|\tau =_\tau (r)\widehat{r},s_z|n,l,j;\mathrm{\Gamma },\nu ,k.`$ (156)
The orientational (spin-orbital) functions are linear combinations of real SAFs and spinors,

$`\widehat{r}|n,l,j;\mathrm{\Gamma },\nu ,k={\displaystyle c_\tau (\lambda ,s_z)S_\lambda (\mathrm{\Omega })\zeta (s_z)}.`$ (157)
In the case of the spherical symmetry $`\lambda =(l,m)`$, $`\tau =(n,j,m_j)`$ and $`c_\tau (\lambda ,s_z)`$ are Clebsch-Gordan (or Wigner) coefficients  but for other point groups these coefficients are not well known. We have derived $`c_\tau (\lambda ,s_z)`$ for the cubic double group $`O_h^{}`$ in Appendix C. Functions $`\widehat{x}|\tau \widehat{r},s_z|\tau `$ are independent of the radial part $`_{n,j,l}(r)`$, which is assumed to be the same for all states belonging to the same $`n`$, $`j`$ and $`l`$. The electronic states $`|n,l,j;\mathrm{\Gamma },\nu ,k`$ distinguished by the index $`k`$ have the same energy, $`ϵ(l,j;\mathrm{\Gamma },\nu )=ϵ(\tau )`$. We next consider the multipole expansion of $`\rho _\tau (\stackrel{}{r})`$, Eq. (154), and find

$`\rho _\tau (\stackrel{}{r})={\displaystyle \underset{\mathrm{\Lambda }}{}}_\tau ^2(r)c_\mathrm{\Lambda }(\tau ,\tau )S_\mathrm{\Lambda }(\widehat{r}).`$ (158)
Since the core density is invariant of the point group, only the full symmetrical irreducible representations has to be considered, with $`\mathrm{\Gamma }_\mathrm{\Lambda }=A_{1g}`$. The coefficients $`c_\mathrm{\Lambda }(\tau ,\tau ^{})`$ are given by

$`c_\mathrm{\Lambda }(\tau ,\tau ^{})={\displaystyle \underset{s_z=1/2}{\overset{+1/2}{}}}{\displaystyle \mathrm{\Omega },s_z|\tau ^{}S_\mathrm{\Lambda }(\mathrm{\Omega })\mathrm{\Omega },s_z|\tau ^{}𝑑\mathrm{\Omega }}.`$ (159)
(We recall that $`\mathrm{\Omega }`$ or $`\widehat{r}`$ stand for the two polar angles $`(\mathrm{\Theta },\varphi )`$.)

A2. The density of the valence electrons reads

$`\rho _{val}(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{k},\alpha }{}}n_{\stackrel{}{k}\alpha }|x|\stackrel{}{k}\alpha |^2,`$ (160)
where $`x|\stackrel{}{k}\alpha `$ is the wave function of a delocalized electron with the wave vector $`\stackrel{}{k}`$ and the band index $`\alpha `$; $`n_{\stackrel{}{k}\alpha }`$ is the occupation number,

We expand the electron density of itinerant electrons in multipole series inside the spheres. In the following we will use the LAPW basis functions, Appendix B. Using Eqs. (B11) we find that the local density at a site $`\stackrel{}{n}`$ is given by

$`\rho _{val}(\stackrel{}{r}(\stackrel{}{n}))={\displaystyle \underset{\mathrm{\Lambda }}{}}\rho _\mathrm{\Lambda }(r(\stackrel{}{n}))S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n})),`$ (161)
where

$`\rho _\mathrm{\Lambda }(r)={\displaystyle \underset{l_1l_2}{}}u_{l_1}^p(r)u_{l_2}^t(r)\rho _\mathrm{\Lambda }^{(l_1,p)(l_2,t)}.`$ (162)
Here and below summation is understood over the repeated indices $`p`$ and $`t`$, and

$`\rho _\mathrm{\Lambda }^{(l_1,p)(l_2,t)}={\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k},\alpha }{}}n_{\stackrel{}{k},\alpha }c_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}(\stackrel{}{k},\alpha ;\stackrel{}{k},\alpha ).`$ (163)
(See Appendix B for definitions.) The coefficients $`c_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}(\stackrel{}{k},\alpha ;\stackrel{}{k},\alpha )`$ are given by

$`c_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}(\stackrel{}{k},\alpha ;\stackrel{}{k},\alpha )={\displaystyle \underset{\lambda _1(l_1)}{}}{\displaystyle \underset{\lambda _2(l_2)}{}}{\displaystyle \underset{s_z}{}}`$ (164)
$`\times [\gamma A]_{\lambda _1}^p(\stackrel{}{k},\alpha ,s_z)[\gamma A]_{\lambda _2}^t(\stackrel{}{k},\alpha ,s_z)c_\mathrm{\Lambda }(\lambda _1,\lambda _2).`$ (165)
We recall that $`\lambda =(l,\mathrm{\Gamma },\nu ,k)`$ (Appendix A) and the summations in (165) are performed over the subindices $`\mathrm{\Gamma },\nu ,k`$ of $`\lambda (l_1)`$ and $`\lambda (l_2)`$ within the manifolds $`l_1`$ and $`l_2`$. Finally,

$`c_\mathrm{\Lambda }(\lambda _1,\lambda _2)={\displaystyle S_{\lambda _1}(\mathrm{\Omega })S_\mathrm{\Lambda }(\mathrm{\Omega })S_{\lambda _2}(\mathrm{\Omega })𝑑\mathrm{\Omega }}.`$ (166)
These coefficients can be tabulated before the self-consistent-field HFR procedure. The density of conduction electrons stays invariant under all symmetry operations, which means that in expansion (166) we consider only irreducible representations of $`A_{1g}`$ symmetry, $`\mathrm{\Gamma }_\mathrm{\Lambda }=A_{1g}`$. In such case the nonzero coefficients in (166) can occur only if (1) $`\mathrm{\Gamma }_{\lambda _1}=\mathrm{\Gamma }_{\lambda _2}`$ or (2) both $`\mathrm{\Gamma }_{\lambda _1}`$ and $`\mathrm{\Gamma }_{\lambda _2}`$ belong to the $`A_{1g}`$ irreducible representation.

In the interstitial region $`\rho _{val}(\stackrel{}{R})`$ is expanded in Fourier series,

$`\rho _{val}(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{K}}{}}\rho _{val}(\stackrel{}{K})e^{i\stackrel{}{K}\stackrel{}{R}},`$ (167)
where from Eqs. (B1) and (B9) we find

$`\rho _{val}(\stackrel{}{K})`$ $`=`$ $`{\displaystyle \frac{1}{Nv}}{\displaystyle \underset{\stackrel{}{k},\alpha }{}}{\displaystyle \underset{\stackrel{}{K}^{}}{}}{\displaystyle \underset{s_z}{}}n_{\stackrel{}{k},\alpha }`$ (169)
$`\times \gamma _\stackrel{}{K}^{}^{}(\stackrel{}{k},\alpha ,s_z)\gamma _{\stackrel{}{K}+\stackrel{}{K}^{}}(\stackrel{}{k},\alpha ,s_z).`$
### B Direct Coulomb matrix elements of core states

If one wants to use the Hartree-Fock-Roothaan method then the problem of the radial dependence of core states arises. In LDA it is solved quite naturally since it is possible to introduce a self-consistent spherically symmetric potential which includes an average exchange term. Nonspherical contributions are small and usually omitted. Then the radial components are obtained through the solution of the Schrödinger (or Dirac) equation in the potential. However, there is no such convenience in the Hartree-Fock approach since the exchange contribution generally can not be reduced to an effective single particle potential. In the HF method one routinely uses radial dependencies of Gaussian or Slater-type (GT or ST). An alternative choice of complete radial basis functions is given in Appendix D. In all cases the radial part is approximated as

$`_\tau (r)={\displaystyle \underset{\eta }{}}\gamma _{\tau ,\eta }u_\eta ^\tau (r),`$ (170)
where $`u_\eta `$ is a radial basis function (GT, ST or other), while $`\eta `$ stands for the parameters specifying the function.

The orientational vectors $`\widehat{x}|j,l,m_j`$ are usually spherical harmonic spinors.  In a crystal, the spherical symmetry is reduced, Eq. (155), and we replace the spinors by their symmetry adapted combinations $`\widehat{x}|\tau `$ as we discussed in Sec. IV A1 and Appendix C. Then the basis function of a core state $`\tau `$ reads

$`\chi _\eta ^\tau (x)x|\tau ,\eta =u_\eta ^\tau (r)\widehat{x}|\tau .`$ (171)
Since the orientational vector $`\widehat{x}|\tau `$ is independent of $`\eta `$, the basis functions are distinguished only by the radial component $`u_\eta ^\tau (r)`$.

In order to calculate the matrix elements of the direct Coulomb interaction for each atomic orbital $`\tau `$ we introduce the Hartree potential so that no self-interaction occurs,

$`V^\tau (\stackrel{}{r})=V(\stackrel{}{R})𝒱_{}^{}{}_{}{}^{\tau }(\stackrel{}{r})=V_0^\tau (\stackrel{}{r})+V_M^\tau (\stackrel{}{r}),`$ (172)
where $`V_0^\tau `$ and $`V_M^\tau `$ refer to the spherical and the nonspherical components of $`V^\tau `$, respectively. First we consider the spherically symmetric component,

$`V_0^\tau (\stackrel{}{r})=V_0(\stackrel{}{R})𝒱_{}^{}{}_{0}{}^{\tau }(\stackrel{}{r}).`$ (173)
We get

$`\tau ,\eta |_0^{Coul}|\tau ,\eta ^{}`$ $`=`$ $`\tau ,\eta |T|\tau ,\eta ^{}`$ (174)
$`+`$ $`{\displaystyle _0^R}u_\eta ^\tau (r)u_\eta ^{}^\tau (r)V_0^\tau (r)r^2𝑑r.`$ (175)
Here $`\tau ,\eta |T|\tau ^{},\eta ^{}`$ are the kinetic energy integrals which are well known for GT (ST) orbitals. For the alternative set of basis radial functions (Appendix D) instead of (175) the matrix elements are given by Eq. (D3).

In order to compute the matrix elements of $`V_M^\tau `$ we partition $`V(\stackrel{}{R})`$ in two parts, Eq. (54). For $`U_S`$ (the potential of the charges inside the sphere $`S(\stackrel{}{n})`$) we distinguish further two contributions, from the itinerant and the core electrons, Eq. (149). Using equations derived in Sec. II, we obtain

$`\tau ,\eta |_M^{Coul}|\tau ,\eta ^{}`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}c_\mathrm{\Lambda }(\tau ,\tau )(D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(out)`$ (176)
$`+`$ $`D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(val)+D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(core)),`$ (177)
where

$`D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(out)={\displaystyle \frac{V_\mathrm{\Lambda }^{out}}{R^l}}q_l^{\eta \eta ^{}}.`$ (178)
Here $`V_\mathrm{\Lambda }^{out}`$ is given by (60) and the multipolar moments $`q_l^{\eta \eta ^{}}`$ are

$`q_l^{\eta \eta ^{}}={\displaystyle _0^R}u_\eta ^\tau (r)u_\eta ^{}^\tau (r)r^{l+2}𝑑r.`$ (179)
For the contribution from the extended states, $`D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(val)`$, we obtain

$`D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(val)={\displaystyle \underset{l_1l_2}{}}C_l(u_{l_1}^pu_{l_2}^t|u_\eta ^\tau u_\eta ^{}^\tau )\rho _\mathrm{\Lambda }^{(l_1,p)(l_2,t)},`$ (180)
where the two-fold radial integrals $`C_l(\mathrm{})`$ are defined by (E2) and $`\rho _\mathrm{\Lambda }^{(l_1,p)(l_2,t)}`$ is defined by Eqs. (163), (165). The contribution from the other core electrons reads

$`D_\mathrm{\Lambda }^{\eta ,\eta ^{}}(core)={\displaystyle \underset{\tau ^{}}{}}C_l(u_\eta ^\tau u_\eta ^{}^\tau |_\tau ^{}_\tau ^{})c_\mathrm{\Lambda }(\tau ^{},\tau ^{})n_\tau ^{}^\tau .`$ (181)
The last sum runs over all occupied states except $`\tau `$ (no self-interaction), i.e. $`n_\tau ^\tau =0`$ and $`n_\tau ^\tau ^{}=1`$ if the core state $`\tau ^{}\tau `$ is occupied and zero otherwise.

In cubic crystals only $`f`$ and $`d`$ shells are split by CEF effects, Appendix C. The CEF interaction of localized electrons with the delocalized ones (second part of (177)) has been discussed in a number of papers, Refs.  and calculated in Refs. . The important result here is that we have obtained all CEF interactions. In particular, besides those considered in Refs. , we include the CEF effects from the rest of the crystal  (first contribution, Eq. (177)) and from the other core shells with nonspherical density (third part of (177)).

### C Direct Coulomb matrix elements of extended states

Here we consider matrix elements of $`^{Coul}`$, Eq. (137), for a conduction electron $`|d|\stackrel{}{k},\alpha `$. Usually in the LAPW method the basis functions $`\chi _{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{R})=\stackrel{}{R}|\stackrel{}{k},\stackrel{}{K}`$ (see Appendix B) are defined in an effective potential $`V_0^{LDA}=V_0+V_0^{exc}`$ which includes the direct Coulomb interaction $`V_0`$ and the LDA exchange potential $`V_0^{exc}`$.  The potential $`V_0^{LDA}`$ is spherically symmetric inside the “muffin-tin” (MT) spheres and is a constant in the interstitial region. Notice that $`V_0^{LDA}`$ is used only for the construction of the basis functions $`\chi _{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{R})`$, Eq. (B1). Having defined $`\chi _{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{R})`$, in principle one can calculate the matrix elements of a general potential $`V^{tot}(\stackrel{}{R})=V^{Coul}(\stackrel{}{R})+V^{exc}(\stackrel{}{R})`$ (the procedure is known as FP-LAPW method ). Below we follow the same approach, but instead of $`V_0^{LDA}`$ for construction of the basis functions we consider only the electrostatic potential $`V_0^d(\stackrel{}{R})`$ without any exchange,

$`V_0^d(\stackrel{}{R})=V_0(\stackrel{}{R})𝒱_{}^{}{}_{0}{}^{d}(\stackrel{}{R}).`$ (182)
Here $`𝒱_{}^{}{}_{0}{}^{d}(\stackrel{}{R})`$ is the spherically symmetric Coulomb potential due to the electron $`d`$. As we will see later in section V we can omit $`𝒱_{}^{}{}_{}{}^{d}(\stackrel{}{R})`$ for any conduction electron, because the corresponding matrix elements, Eq. (279), decreases as

$`\stackrel{}{k},\stackrel{}{K}^{}|𝒱_{}^{}{}_{}{}^{d}|\stackrel{}{k},\stackrel{}{K}{\displaystyle \frac{1}{N}},`$ (183)
and vanish in the limit $`N\mathrm{}`$. Therefore, constructing the LAPW basis functions $`\chi _{\stackrel{}{k},\stackrel{}{K}}`$, we can use the potential $`V_0(\stackrel{}{R})`$, i.e. $`V_0^d(\stackrel{}{R})=V_0(\stackrel{}{R})`$. Next step is to calculate the direct matrix elements for the conduction electron $`d`$ in the full Coulomb potential $`V(\stackrel{}{R})`$. Following the method described in Sec. II we write $`V(\stackrel{}{R})=V_0(\stackrel{}{R})+V_M(\stackrel{}{R})`$ for the potential inside the spheres ($`V_0`$ and $`V_M`$ are the spherically symmetric part and the contribution due to the other multipoles, correspondingly). In the interstitial region we use the Fourier expansion of $`V(\stackrel{}{R})`$. Then the Hartree-Fock operator $`^{Coul}`$ is separated into three parts,

$`^{Coul}=_0^{Coul}+_M^{Coul}+_I^{Coul},`$ (184)
where $`_0^{Coul}`$ comprises the kinetic energy and the spherically symmetric potential $`V_0(\stackrel{}{R})`$, $`F_M^{Coul}`$ accounts for the other multipole terms while $`_I^{Coul}`$ stands for the electrostatic interaction in the interstitial region. Starting with $`_0^{Coul}`$ we arrive at the standard expressions for the matrix elements of the LAPW method,

$`\stackrel{}{p},\stackrel{}{K}^{},s_z^{}|_0^{Coul}+K_I|\stackrel{}{k},\stackrel{}{K},s_z=\delta _{\stackrel{}{p},\stackrel{}{k}}\delta _{s_zs_z^{}}H_{\stackrel{}{K}^{},\stackrel{}{K}}(\stackrel{}{k})`$ (185)
where $`H_{\stackrel{}{K}^{},\stackrel{}{K}}(\stackrel{}{k})`$ is given by Eq. (15) of Ref. . (The overlap matrix is given by Eq. (13a) of Ref. .)

The other contributions of the general Coulomb potential $`V^{Coul}(\stackrel{}{R})`$ follow from equations of Sec. II. After some algebra, we obtain:

$`\stackrel{}{p},\stackrel{}{K}^{},s_z^{}|_M^{Coul}+_I^{Coul}|\stackrel{}{k},\stackrel{}{K},s_z=\delta _{\stackrel{}{k}\stackrel{}{p}}\delta _{s_zs_z^{}}`$ (186)
$`\times \left({\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \underset{l_1l_2}{}}B_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}\overline{c}_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}+B_I\right).`$ (187)
Here $`B_I`$ arises due to the interstitial contribution,

$`B_I={\displaystyle \underset{\stackrel{}{P}0}{}^{}}V(\stackrel{}{P})𝒪(\stackrel{}{K}\stackrel{}{K}^{}+\stackrel{}{P}),`$ (188)
where $`V(\stackrel{}{P})`$ stands for the Fourier coefficients, Eq. (49). $`𝒪(\stackrel{}{P}^{})`$ is given by (90) and

$`\overline{c}_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}={\displaystyle \underset{\lambda _1(l_1)}{}}{\displaystyle \underset{\lambda _2(l_2)}{}}A_{\lambda _1}^p(\stackrel{}{k},\stackrel{}{K}^{})A_{\lambda _2}^t(\stackrel{}{k},\stackrel{}{K})c_\mathrm{\Lambda }(\lambda _1,\lambda _2).`$ (189)
(190)
To condense notations we shall use $`L`$ for $`(l,p)`$. In (186)

$`B_\mathrm{\Lambda }^{L_1L_2}=B_\mathrm{\Lambda }^{L_1L_2}(out)+B_\mathrm{\Lambda }^{L_1L_2}(val)+B_\mathrm{\Lambda }^{L_1L_2}(core),`$ (191)
where the index $`out`$ stands for the contribution from the other spheres,

$`B_\mathrm{\Lambda }^{L_1L_2}(out)={\displaystyle \frac{V_\mathrm{\Lambda }^{out}}{R^l}}q_l^{L_1L_2}.`$ (192)
Here $`V_\mathrm{\Lambda }^{out}`$ is given by (60) and the multipolar moments $`q_\mathrm{\Lambda }^{L_1L_2}`$ are

$`q_l^{(l_1,p)(l_2,t)}={\displaystyle _0^R}u_{l_1}^p(r)u_{l_2}^t(r)r^{l+2}𝑑r.`$ (193)
For the contribution from the delocalized electrons we obtain

$`B_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}(val)={\displaystyle \underset{l_1^{}l_2^{}}{}}C_l(u_{l_1}^pu_{l_2}^t|u_{l_1^{}}^p^{}u_{l_2^{}}^t^{})\rho _\mathrm{\Lambda }^{(l_1^{},p^{})(l_2^{},t^{})},`$ (194)
where $`\rho _\mathrm{\Lambda }^{L_1L_2}`$ and $`C_l(\mathrm{})`$ are given by Eqs. (163) and (E2), respectively.

For the contribution from the closed shell core electrons we have

$`B_\mathrm{\Lambda }^{(l_1,p)(l_2,t)}(core)={\displaystyle \underset{\tau }{}}C_l(u_{l_1}^pu_{l_2}^t|_\tau _\tau )c_\mathrm{\Lambda }(\tau ,\tau )n_\tau ,`$ (195)
where again the two-fold radial integral $`C_l`$ is defined by Eq. (E2).

The non-spherical components of the total Coulomb potential in (186) represent crystal field like effects for conduction electrons.

## V Calculations of the Exchange Matrix Elements

In this section we calculate the exchange matrix elements of the Hartree-Fock method. The most general expression for exchange is

$`a|^{exc}|b={\displaystyle \underset{c}{}}{\displaystyle \chi _a^{}(x_1)\psi _c(x_1)\frac{1}{R_{12}}\chi _b(x_2)\psi _c^{}(x_2)𝑑x_1𝑑x_2}.`$ (196)
(197)
Here both indices $`a`$ and $`b`$ refer to basis wave functions of an electronic state $`d`$ (which can be either a conduction or a localized state), Eq. (129); $`\psi _c`$ stand for the estimated wave functions of conduction and localized electrons obtained from a previous iteration of the HFR self-consistent procedure. The summation is understood over all occupied states $`c`$. As before (Sec. III), $`x=(\stackrel{}{R},s_z)`$ and the integration (197) includes summation over two spin components $`s_z`$. We will calculate the matrix element (197) in two steps: 1) we construct an auxiliary Coulomb potential

$`V_{cb}(\stackrel{}{R})={\displaystyle \frac{\rho _{cb}(\stackrel{}{R}^{})}{|\stackrel{}{R}\stackrel{}{R}^{}|}𝑑\stackrel{}{R}^{}},`$ (198)
which corresponds to the “exchange” density

$`\rho _{cb}(\stackrel{}{R})={\displaystyle \underset{s_z=1/2}{\overset{1/2}{}}}\psi _c^{}(\stackrel{}{R},s_z)\chi _b(\stackrel{}{R},s_z).`$ (199)
($`\psi _c`$ has two spin components due to the spin-orbit interaction, Eq. (135).); 2) We calculate the matrix element of exchange as

$`a|^{exc}|b={\displaystyle \underset{c}{}}{\displaystyle V_{cb}(\stackrel{}{R})\rho _{ca}^{}(\stackrel{}{R})𝑑\stackrel{}{R}}.`$ (200)
We want to stress that the “exchange” potential $`V_{cb}`$ and the “exchange” density $`\rho _{cb}`$ are technical quantities here, which are employed only for the calculation of the matrix element of exchange, Eq. (200). They should not be confused with the effective exchange potential which is derived and widely used in the local density approximation.

### A Exchange for localized electrons

Here we calculate the exchange for a localized state $`\tau `$ sited at $`\stackrel{}{n}`$, $`i.e.`$ $`d\tau `$. Index $`c`$ can refer either to a conduction state or to another core state $`\tau ^{}`$. First we consider $`c`$ as a conduction state, $`i.e.`$ $`|c|\stackrel{}{k},\alpha `$. The radial wave function of the localized electron $`\tau `$ is expanded in terms of $`u_\eta ^\tau `$, Eq. (170), so that $`|a`$ and $`|b`$ stand for $`|\tau ,\eta `$. Since the state $`\tau `$ is confined inside $`S(\stackrel{}{n})`$, for the calculation of the exchange (200) we need to know the density $`\rho _{cb}`$ and the potential $`V_{cb}`$, Eqs. (199) and (198), only inside the same sphere. For $`\rho _{cb}`$ we get

$`\rho _{(\stackrel{}{k},\alpha )(\tau ,\eta )}(\stackrel{}{r}(\stackrel{}{n}))={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{\mathrm{\Lambda }}{}}\rho _\mathrm{\Lambda }^{(\stackrel{}{k},\alpha )\tau }(r(\stackrel{}{n}))S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n})),`$ (201)
(202)
where

$`\rho _\mathrm{\Lambda }^{(\stackrel{}{k},\alpha )(\tau ,\eta )}(r)`$ $`=`$ $`{\displaystyle \underset{\lambda ^{}}{}}{\displaystyle \underset{s_z}{}}c_\mathrm{\Lambda }(\tau ,s_z;\lambda ^{})[\gamma A]_\lambda ^{}^p(\stackrel{}{k},\alpha ,s_z)`$ (204)
$`\times u_l^{}^p(r)u_\eta ^\tau (r).`$
The coefficients $`c_\mathrm{\Lambda }(\tau ,s_z;\lambda ^{})`$ are given by

$`c_\mathrm{\Lambda }(\tau ,s_z;\lambda )={\displaystyle \mathrm{\Omega },s_z|\tau S_\mathrm{\Lambda }(\mathrm{\Omega })S_\lambda (\mathrm{\Omega })𝑑\mathrm{\Omega }}.`$ (205)
From (202) the exchange between $`|\tau `$ and $`|\stackrel{}{k},\alpha `$ is found as

$`\tau ,\eta |^{exc}(\stackrel{}{k},\alpha )|\tau ,\eta ^{}={\displaystyle \frac{1}{N}}D^{exc}(\stackrel{}{k},\alpha ),`$ (206)
where

$`D^{exc}(\stackrel{}{k},\alpha )={\displaystyle \underset{\mathrm{\Lambda }}{}}{\displaystyle \underset{\lambda _1\lambda _2}{}}{\displaystyle \underset{s_z,s_z^{}}{}}c_\mathrm{\Lambda }(\tau ,s_z;\lambda _1)c_\mathrm{\Lambda }(\tau ,s_z^{};\lambda _2)`$ (207)
$`\times [\gamma A]_{\lambda _1}^p(\stackrel{}{k},\alpha ,s_z)[\gamma A]_{\lambda _2}^p^{}(\stackrel{}{k},\alpha ,s_z^{})C_l(u_{l_1}^pu_\eta ^\tau |u_{l_2}^p^{}u_\eta ^{}^\tau ).`$ (208)
Here the integrals $`C_l(\mathrm{})`$ are given by Eq. (E2). The exchange with all extended states reads

$`\tau ,\eta |^{exc}(c)|\tau ,\eta ^{}={\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k},\alpha }{}}D^{exc}(\stackrel{}{k},\alpha )n_{\stackrel{}{k},\alpha }.`$ (209)
Next we consider the exchange (200) between $`\tau `$ (as before, $`\tau d`$) and a core electron $`\tau ^{}`$ localized at the same site $`\stackrel{}{n}`$. $`\tau ^{}`$ is now our reference state $`c`$, i.e. $`c=\tau ^{}`$. Proceeding analogously, we find that the multipole expansion of the “exchange” density $`\rho _{cb}`$ is

$`\rho _{\tau ^{}(\tau ,\eta )}(\stackrel{}{r}(\stackrel{}{n}))={\displaystyle \underset{\mathrm{\Lambda }}{}}\rho _\mathrm{\Lambda }^{\tau ^{}(\tau ,\eta )}(r(\stackrel{}{n}))S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n})),`$ (210)
where

$`\rho _\mathrm{\Lambda }^{\tau ^{}(\tau ,\eta )}(r)=c_\mathrm{\Lambda }(\tau ,\tau ^{})_\tau ^{}(r)u_\eta ^\tau (r).`$ (211)
The exchange integral is

$`\tau ,\eta |^{exc}(core)|\tau ,\eta ^{}={\displaystyle \underset{\tau ^{}}{}}n_\tau ^\tau ^{}D^{exc}(\tau ^{}),`$ (213)
where

$`D^{exc}(\tau ^{})={\displaystyle \underset{\mathrm{\Lambda }}{}}|c_\mathrm{\Lambda }(\tau ,\tau ^{})|^2C_l(_\tau ^{}u_\eta ^\tau |_\tau ^{}u_\eta ^{}^\tau ).`$ (214)
### B Exchange between an extended state with the localized states

In this subsection we consider $`d`$ as an extended electron state, i.e. $`d=(\stackrel{}{k},\alpha )`$, while $`c`$ refers to a core state $`\tau `$ localized inside a sphere $`S(\stackrel{}{n})`$, $`c=\tau (\stackrel{}{n})`$. We expand the wave function $`\stackrel{}{R}|\stackrel{}{k},\alpha `$ in terms of $`\stackrel{}{R},s_z|\stackrel{}{k},\stackrel{}{K}=\stackrel{}{R}|\stackrel{}{k},\stackrel{}{K}\zeta (s_z)`$, Appendix B. Therefore, $`a,b=(\stackrel{}{k},\stackrel{}{K},s_z)`$. The exchange density $`\rho _{cb}`$, Eq. (199), is located inside the sphere $`S(\stackrel{}{n})`$,

$`\rho _{\tau (\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{r})={\displaystyle \frac{1}{\sqrt{N}}}e^{i\stackrel{}{k}\stackrel{}{X}(\stackrel{}{n})}{\displaystyle \underset{\mathrm{\Lambda }}{}}{\displaystyle \underset{\lambda ^{}}{}}(_\tau (r)u_l^{}^p(r)`$ (215)
$`\times c_\mathrm{\Lambda }(\tau ,s_z;\lambda ^{})A_\lambda ^{}^p(\stackrel{}{k},\stackrel{}{K}))S_\mathrm{\Lambda }(\widehat{r}),`$ (216)
where the coefficients $`c_\mathrm{\Lambda }(\tau ,s_z;\lambda ^{})`$ are given by Eq. (205). Again, using the multipole expansion (216) we calculate the exchange (200) and obtain

$`\stackrel{}{k},\stackrel{}{K},s_z|^{exc}(\tau (\stackrel{}{n}))|\stackrel{}{k},\stackrel{}{K}^{},s_z^{}={\displaystyle \frac{1}{N}}B^{exc}(\tau ),`$ (217)
where

$`B^{exc}(\tau )={\displaystyle \underset{\mathrm{\Lambda }}{}}{\displaystyle \underset{\lambda _1\lambda _2}{}}c_\mathrm{\Lambda }(\tau ,s_z;\lambda _1)c_\mathrm{\Lambda }(\tau ,s_z^{};\lambda _2)`$ (218)
$`\times A_{\lambda _1}^p(\stackrel{}{k},\stackrel{}{K})A_{\lambda _2}^t(\stackrel{}{k},\stackrel{}{K}^{})C_l(_\tau u_{l_1}^p|_\tau u_{l_2}^t).`$ (219)
Notice that the latter result is independent of $`\stackrel{}{n}`$ and the matrix elements (217) can be nonzero even for off-diagonal spin functions ($`s_z=1/2`$, $`s_z^{}=1/2`$ or vice versa) due to the spin-orbit interaction. In order to obtain the exchange with all localized electrons, we sum (217) over all sites $`\stackrel{}{n}`$ and all occupied core electron states $`|\tau `$, and find that

$`\stackrel{}{k},\stackrel{}{K},s_z|^{exc}(core)|\stackrel{}{k},\stackrel{}{K}^{},s_z^{}={\displaystyle \underset{\tau }{}}B^{exc}(\tau ).`$ (220)
### C Exchange between extended states

The calculation of exchange between two delocalized electrons is much more involved and quite laborious. As the state $`c`$ in Eqs. (197)-(200) we consider now a conduction state $`|\stackrel{}{p},\beta `$ with the wave function $`\psi _c(\stackrel{}{R})=x|\stackrel{}{p},\beta `$. As before, $`d=(\stackrel{}{k},\alpha )`$ and the wave function of the $`d`$th electron, $`\stackrel{}{R}|\stackrel{}{k},\alpha `$, is expanded in terms of LAPW basis functions $`\stackrel{}{R},s_z|\stackrel{}{k},\stackrel{}{K}=\stackrel{}{R}|\stackrel{}{k},\stackrel{}{K}\zeta (s_z)`$ (Appendix B) labeled by indices $`a,b=(\stackrel{}{k},\stackrel{}{K},s_z)`$. In this subsection we first calculate the matrix elements of $`d`$th electron exchanged with the extended state $`c`$. Then by summing over all occupied extended states $`|c|\stackrel{}{p},\beta `$ we will be able to compute the matrix element of exchange between $`d`$ and the other conduction electrons.

Inside a sphere $`S(\stackrel{}{n})`$ the “exchange” density $`\rho _{cb}`$, Eq. (199), is given by

$`\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{r})={\displaystyle \frac{1}{N}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}{\displaystyle \underset{\mathrm{\Lambda }}{}}\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r(\stackrel{}{n}))S_\mathrm{\Lambda }(\widehat{r}(\stackrel{}{n})),`$ (221)
(222)
where $`\stackrel{}{q}=\stackrel{}{p}\stackrel{}{k}`$ and

$`\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)`$ $`=`$ $`{\displaystyle \underset{l_1l_2}{}}{\displaystyle \underset{\lambda _1\lambda _2}{}}c_\mathrm{\Lambda }(\lambda _1,\lambda _2)[\gamma A]_{\lambda _1}^p(\stackrel{}{p},\beta ,s_z)`$ (224)
$`\times A_{\lambda _2}^t(\stackrel{}{k},\stackrel{}{K})u_{l_1}^p(r)u_{l_2}^t(r).`$
Outside the spheres $`\rho _{cb}`$ is expanded in plane waves:

$`\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})={\displaystyle \frac{1}{Nv}}e^{i\stackrel{}{q}\stackrel{}{R}}{\displaystyle \underset{\stackrel{}{K}^{}}{}}\gamma _{\stackrel{}{K}^{}+\stackrel{}{K}}^{}(\stackrel{}{p},\beta ,s_z)e^{i\stackrel{}{K}^{}\stackrel{}{R}}.`$ (225)
Proceeding as in Sec. II we continue the plane wave representation (225) inside the spheres and then subtract it out from Eq. (222). This procedure renormalizes the multipole radial functions $`\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}`$ of (224), which now are given by

$`\rho _{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)=\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r){\displaystyle \frac{4\pi i^l}{v}}`$ (226)
$`\times {\displaystyle \underset{\stackrel{}{K}^{}}{}}j_l(|\stackrel{}{K}^{}+\stackrel{}{q}|r)S_\mathrm{\Lambda }(\stackrel{}{K}^{}+\stackrel{}{q})\gamma _{\stackrel{}{K}^{}+\stackrel{}{K}}^{}(\stackrel{}{p},\beta ,s_z).`$ (227)
Here $`S_\mathrm{\Lambda }(\stackrel{}{K}^{}+\stackrel{}{q})`$ is understood as the symmetry adapted function $`S_\mathrm{\Lambda }(\mathrm{\Omega }_{\stackrel{}{K}^{}+\stackrel{}{q}})`$ of the polar angles $`\mathrm{\Omega }_{\stackrel{}{K}^{}+\stackrel{}{q}}`$ defined by the vector $`\stackrel{}{K}^{}+\stackrel{}{q}`$.

Notice that $`\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}`$ is not a periodic function of $`\stackrel{}{R}`$. A translation by $`\stackrel{}{t}=\stackrel{}{X}(\stackrel{}{n})`$ transforms it as

$`\stackrel{}{t}\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}).`$ (229)
This transformation is not identical for $`\stackrel{}{q}0`$. From Eq. (198) we find that the same transformational law holds for the corresponding potential, i.e.

$`\stackrel{}{t}V_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}V_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}).`$ (230)
As a consequence of Bloch’s theorem we obtain

$`\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=e^{i\stackrel{}{q}\stackrel{}{R}}\stackrel{~}{\rho }_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}),`$ (232)
$`V_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=e^{i\stackrel{}{q}\stackrel{}{R}}\stackrel{~}{V}_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}),`$ (233)
where $`\stackrel{~}{\rho }_{cb}`$ and $`\stackrel{~}{V}_{cb}`$ are periodic functions, i.e.

$`\stackrel{}{t}\stackrel{~}{\rho }_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=\stackrel{~}{\rho }_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}),`$ (235)
$`\stackrel{}{t}\stackrel{~}{V}_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R})=\stackrel{~}{V}_{c(\stackrel{}{k},\stackrel{}{K},s_z)}(\stackrel{}{R}).`$ (236)
Therefore, in the interstitial region $`\stackrel{~}{\rho }_{cb}`$ and $`\stackrel{~}{V}_{cb}`$ are expanded in Fourier series, while the initial functions $`\rho _{cb}`$ and $`V_{cb}`$ are expanded in terms of $`\mathrm{exp}(i(\stackrel{}{q}+\stackrel{}{K}^{})\stackrel{}{R})`$, where $`\stackrel{}{K}^{}`$ is a reciprocal lattice vector. The factor $`\mathrm{exp}(i\stackrel{}{q}\stackrel{}{R})`$ or, more precisely, $`\mathrm{exp}(i\stackrel{}{q}\stackrel{}{X})`$, will be also present for solutions inside the spheres. Indeed, if one knows a solution inside a sphere, say $`S(\stackrel{}{n}=0)`$ at the origin, then by means of Eq. (230) it is easy to generate the solution inside any other sphere. Notice, however, that the factor $`\mathrm{exp}(i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n}))`$ cancels in the final expression (200). As a result, it is easily to generalize the method of Sec. II for the present consideration.

First we consider the potential in the interstitial region,

$`V_e(\stackrel{}{R})={\displaystyle \frac{1}{N}}e^{i\stackrel{}{q}\stackrel{}{R}}{\displaystyle \underset{\stackrel{}{K}^{}}{}}V_e(\stackrel{}{K}^{})e^{i\stackrel{}{K}^{}\stackrel{}{R}}.`$ (237)
As in Sec. II we distinguish three contributions there,

$`V_e(\stackrel{}{K}^{})=V_e^I(\stackrel{}{K}^{})+V_e^{S,0}(\stackrel{}{K}^{})+V_e^{S,M}(\stackrel{}{K}^{}),`$ (238)
where $`V_e^{S,0}(\stackrel{}{K}^{})`$ and $`V_e^{S,M}(\stackrel{}{K}^{})`$ are the Fourier components of the potentials of the monopole and the other higher multipoles of the sphere, respectively. $`V_e^I(\stackrel{}{K}^{})`$ represents a component from the the plane wave expansion (225), i.e.

$`V_e^I(\stackrel{}{K}^{})={\displaystyle \frac{1}{v}}{\displaystyle \frac{4\pi }{|\stackrel{}{K}^{}+\stackrel{}{q}|^2}}\gamma _{\stackrel{}{K}^{}+\stackrel{}{K}}^{}(\stackrel{}{p},\beta ,s_z).`$ (239)
The Fourier components $`V_e^{S,0}(\stackrel{}{K}^{})`$ and $`V_e^{S,M}(\stackrel{}{K}^{})`$ are found as

$`V_e^{S,g}(\stackrel{}{K}^{})={\displaystyle \frac{1}{v}}{\displaystyle _{cell}}e^{i(\stackrel{}{q}+\stackrel{}{K}^{})\stackrel{}{R}}V_{S_g}(\stackrel{}{R})d^3R,`$ (240)
where $`g=0`$, $`M`$, and

$`V_{S_g}(\stackrel{}{R})={\displaystyle \underset{\stackrel{}{n}}{}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}𝒰_{S_g}(\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})).`$ (241)
Here $`𝒰_{S_0}`$ and $`𝒰_{S_M}`$ are the potentials due to the monopole ($`l=0`$) and the other multipoles ($`l1`$) of a single sphere (see Eq. (18)). For the potential $`𝒰_S`$, Eq. (18), we use the equivalent charge distribution given by Eq. (227). The integral (240) is taken over a unit cell of the Bravais lattice. One can show that $`V_e^{S,g}(\stackrel{}{K}^{})`$ is also expressed as

$`V_e^{S,g}(\stackrel{}{K}^{})={\displaystyle \frac{1}{v}}{\displaystyle _{crystal}}e^{i(\stackrel{}{q}+\stackrel{}{K}^{})\stackrel{}{R}}𝒰_{S_g}(\stackrel{}{R})d^3R,`$ (242)
where the integration spans the whole crystal. The representation (242) is then used to compute the Ewald expansion coefficients $`V_e^{S,0}(\stackrel{}{K}^{})`$ and to find $`V_e^{S,M}(\stackrel{}{K}^{})`$ by direct integration. The procedure is the same as in Sec. II.

In order to proceed with the Ewald expansion, we introduce an effective point “exchange charge” inside the spheres,

$`Q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}=\sqrt{4\pi }{\displaystyle _0^R}\rho _{}^{}{}_{0}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)r^2.`$ (243)
Using the orthogonality of the radial functions $`u_l^1`$ and $`u_l^2`$ and the relation $`c_0(\lambda ,\lambda ^{})=\delta _{\lambda \lambda ^{}}/\sqrt{4\pi }`$, we obtain

$`Q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}={\displaystyle \underset{l}{}}{\displaystyle \underset{\lambda }{}}[\gamma A]_\lambda ^p(\stackrel{}{p},\beta ,s_z)A_\lambda ^p(\stackrel{}{k},\stackrel{}{K})N_l^p`$ (244)
$`{\displaystyle \frac{4\pi R^2}{v}}{\displaystyle \underset{\stackrel{}{K}^{}}{}}{\displaystyle \frac{j_1(|\stackrel{}{K}^{}+\stackrel{}{q}|R)}{|\stackrel{}{K}^{}+\stackrel{}{q}|}}\gamma _{\stackrel{}{K}^{}+\stackrel{}{K}}^{}(\stackrel{}{p},\beta ,s_z),`$ (245)
where $`N_l^1=1`$ (normalization of the radial functions of the LAPW method) and

$`N_l^2`$ $`=`$ $`{\displaystyle _0^R}[u_l^2(r)]^2r^2𝑑r.`$ (246)
Although in general $`Q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(R)0`$ in Appendix F we prove that for the present particular case (one atom in the unit cell) this “charge” has to be considered only for $`\stackrel{}{q}=0`$, i.e. when $`\stackrel{}{p}=\stackrel{}{q}`$ but $`\alpha \beta `$. However, the final results will not be affected if one assumes that $`Q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(R)0`$ for $`\stackrel{}{q}0`$. If there are two or more atoms in the unit cell, the “exchange charges” of nonequivalent atoms must be introduced and taken into account even for the case $`\stackrel{}{p}\stackrel{}{q}`$, Appendix F. We shall now proceed further as for the general case. The Ewald expansion coefficients are given by

$`V_e^{S,0}(\stackrel{}{K}^{})={\displaystyle \frac{\pi }{v}}{\displaystyle \frac{4}{|\stackrel{}{K}^{}+\stackrel{}{q}|^2}}e^{|\stackrel{}{K}^{}+\stackrel{}{q}|^2/4G^2}Q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}.`$ (247)
For the higher multipoles we introduce two functions,

$`Q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)={\displaystyle _0^r}\rho _{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r_{}^{}{}_{}{}^{l+2}𝑑r^{},`$ (249)
$`Q_{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)={\displaystyle _r^R}\rho _{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}1}l}𝑑r^{}.`$ (250)
Proceeding then as in section II we find:

$`V_e^{S,M}(\stackrel{}{K}^{})={\displaystyle \frac{(4\pi )^2}{v}}{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}{\displaystyle \frac{(i)^l}{2l+1}}S_\mathrm{\Lambda }(\stackrel{}{K}^{}+\stackrel{}{q})`$ (251)
$`\times \left(A_\mathrm{\Lambda }(|\stackrel{}{K}^{}+\stackrel{}{q}|)+Q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(R){\displaystyle \frac{j_{l1}(|\stackrel{}{K}^{}+\stackrel{}{q}|R)}{|\stackrel{}{K}^{}+\stackrel{}{q}|R^{l1}}}\right),`$ (252)
where

$`A_\mathrm{\Lambda }(|\stackrel{}{K}^{}+\stackrel{}{q}|)`$ $`=`$ $`{\displaystyle _0^R}\left({\displaystyle \frac{Q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)}{r^{l+1}}}+r^lQ_{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)\right)`$ (254)
$`\times j_l(|\stackrel{}{K}^{}+\stackrel{}{q}|r)r^2dr.`$
Eqs. (239), (247) and (252) fully determine the plane wave components $`V_e(\stackrel{}{K}^{})`$, Eq. (238). Thus, we have obtained the potential $`V_e(\stackrel{}{R})`$, Eq. (237), in the interstitial region.

The potential inside a sphere $`S(\stackrel{}{n})`$ is found as

$`V_e(\stackrel{}{r})={\displaystyle \frac{1}{N}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}[V_{e,0}(r)+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}V_{e,\mathrm{\Lambda }}(r)S_\mathrm{\Lambda }(\widehat{r})],`$ (255)
where $`V_{e,0}(r)`$ and $`V_{e,\mathrm{\Lambda }}(r)`$ are smooth functions of $`r`$. As before, we introduce the constants $`V_{e,0}V_{e,0}(R)`$ and $`V_{e,\mathrm{\Lambda }}V_{e,\mathrm{\Lambda }}(R)`$. They are found from the boundary-value problem for the surface of the sphere:

$`V_{e,\mathrm{\hspace{0.17em}0}}={\displaystyle \underset{\stackrel{}{K}^{}}{}}j_0(|\stackrel{}{K}^{}+\stackrel{}{q}|R)V_e(\stackrel{}{K}^{}),`$ (257)
$`V_{e,\mathrm{\Lambda }}=4\pi i^l{\displaystyle \underset{\stackrel{}{K}^{}}{}}j_l(|\stackrel{}{K}^{}+\stackrel{}{q}|R)S_\mathrm{\Lambda }(\stackrel{}{K}^{}+\stackrel{}{q})V_e(\stackrel{}{K}^{}).`$ (258)
(259)
Here $`V_{e,\mathrm{\hspace{0.17em}0}}`$ and $`V_{e,\mathrm{\Lambda }}`$ depend on $`\stackrel{}{q}`$ and have to be recalculated for each extended state $`\stackrel{}{p}=\stackrel{}{k}+\stackrel{}{q}`$. We recall that $`V_e(\stackrel{}{K}^{})`$ is given by the sum (238). Following Sec. II.C we rewrite the potential $`V_e(\stackrel{}{r})`$, Eq. (255), as

$`V_e(\stackrel{}{r}(\stackrel{}{n}))={\displaystyle \frac{1}{N}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}\left[U_{e,S}(\stackrel{}{r}(\stackrel{}{n}))+V_{e,S}^{out}(\stackrel{}{r}(\stackrel{}{n}))\right],`$ (260)
where $`U_{e,S}(\stackrel{}{r}(\stackrel{}{n}))`$ is the potential of the charges located inside the sphere $`S(\stackrel{}{n})`$ and $`V_{e,S}^{out}(\stackrel{}{r})`$ accounts for the potential of the rest of the crystal. The single sphere potential $`U_{e,S}(\stackrel{}{r})`$ is given by Eq. (56) where $`q_0(r)`$, $`q_0^{}(r)`$ are replaced by

$`q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)=\sqrt{4\pi }{\displaystyle _0^r}\rho _0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r_{}^{}{}_{}{}^{2}𝑑r^{},`$ (262)
$`q_{}^{}{}_{0}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)=\sqrt{4\pi }{\displaystyle _r^R}\rho _0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r^{}𝑑r^{},`$ (263)
and $`q_\mathrm{\Lambda }(r)`$, $`q_\mathrm{\Lambda }^{}(r)`$ by

$`q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)={\displaystyle _0^r}\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r_{}^{}{}_{}{}^{l+2}𝑑r^{},`$ (264)
$`q_{}^{}{}_{\mathrm{\Lambda }}{}^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r)={\displaystyle _r^R}\rho _\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(r^{})r_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}1}l}𝑑r^{}.`$ (265)
The potential of the “exchange charges” outside the sphere $`S(\stackrel{}{n})`$ is found by subtracting from Eq. (257,b) the potential of the “charges” inside. Thus, we introduce the constants

$`V_{e,\mathrm{\hspace{0.17em}0}}^{out}=V_{e,\mathrm{\hspace{0.17em}0}}{\displaystyle \frac{q_0^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(R)}{R}},`$ (267)
$`V_{e,\mathrm{\Lambda }}^{out}=V_{e,\mathrm{\Lambda }}{\displaystyle \frac{4\pi }{2l+1}}{\displaystyle \frac{q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{K},s_z)}(R)}{R^{l+1}}}.`$ (268)
Then the potential inside the sphere due to the “exchange charges” outside is

$`V_{e,S}^{out}(\stackrel{}{r})={\displaystyle \frac{1}{N}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}[V_{e,\mathrm{\hspace{0.17em}0}}^{out}+{\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}V_{e,\mathrm{\Lambda }}^{out}\left({\displaystyle \frac{r}{R}}\right)^lS_\mathrm{\Lambda }(\widehat{r})].`$ (269)
(270)
Therefore, Eqs. (260), (270) and (56) where $`q_0(r)`$, $`q_0^{}(r)`$, $`q_\mathrm{\Lambda }(r)`$, $`q_\mathrm{\Lambda }^{}(r)`$ are given by (262-d) fully determine the effective “exchange” potential inside any sphere.

In order to obtain the matrix element of exchange we integrate the “exchange” potential with the “exchange” density $`\rho _{c(\stackrel{}{k},\stackrel{}{K},s_z)}^{}(\stackrel{}{r})`$, Eqs. (222)-(225), and recall that $`c`$ stands for $`(\stackrel{}{p},\beta )`$. From the potentials inside the spheres and in interstices we distinguish three contributions to the exchange,

$`\stackrel{}{k},\stackrel{}{P},s_z^{}|^{exc}(\stackrel{}{p},\beta )|\stackrel{}{k},\stackrel{}{K},s_z={\displaystyle \frac{1}{N}}(B_{out}^{exc}(\stackrel{}{p},\beta )`$ (271)
$`+B_I^{exc}(\stackrel{}{p},\beta )+B_S^{exc}(\stackrel{}{p},\beta )).`$ (272)
Here $`B_{out}^{exc}`$, $`B_S^{exc}`$ and $`B_I^{exc}`$ are the contributions from integrations with $`V_{e,S}^{out}(\stackrel{}{r})`$, $`U_{e,S}(\stackrel{}{r})`$ inside a sphere, and with $`V_e(\stackrel{}{R})`$ in the interstices, respectively. First we carry out the integrations inside a sphere $`S(\stackrel{}{n})`$, and then perform summation over the $`N`$ spheres. As a result we get

$`B_{out}^{exc}(\stackrel{}{p},\beta )={\displaystyle \underset{\mathrm{\Lambda }0}{}^{}}V_{e,\mathrm{\Lambda }}^{out}{\displaystyle \frac{q_\mathrm{\Lambda }^{c(\stackrel{}{k},\stackrel{}{P},s_z^{})}(R)}{R^l}}.`$ (273)
The contribution from the interstitial region gives

$`B_I^{exc}(\stackrel{}{p},\beta )={\displaystyle \underset{\stackrel{}{K}^{},\stackrel{}{P}}{}}V_e(\stackrel{}{K}^{})\gamma _{\stackrel{}{P}+\stackrel{}{K}}(\stackrel{}{p},\beta ,s_z^{})𝒪(\stackrel{}{K}^{}\stackrel{}{P}),`$ (274)
where $`𝒪(\stackrel{}{K}^{\prime \prime })`$ is the overlap integral (90). Finally, for $`B_S^{exc}`$ we get

$`B_S^{exc}(\stackrel{}{p},\beta )={\displaystyle \underset{\mathrm{\Lambda }}{}}{\displaystyle \underset{l_1l_2}{}}{\displaystyle \underset{l_1^{}l_2^{}}{}}{\displaystyle \underset{\lambda _1\lambda _2}{}}{\displaystyle \underset{\lambda _1^{}\lambda _2^{}}{}}c_\mathrm{\Lambda }(\lambda _1,\lambda _2)c_\mathrm{\Lambda }(\lambda _1^{},\lambda _2^{})`$ (275)
$`\times [\gamma A]_{\lambda _1}^p(\stackrel{}{p},\beta ,s_z^{})[\gamma A]_{\lambda _1^{}}^p^{}(\stackrel{}{p},\beta ,s_z)A_{\lambda _2}^t(\stackrel{}{k},\stackrel{}{K})A_{\lambda _2^{}}^t^{}(\stackrel{}{k},\stackrel{}{P})`$ (276)
$`\times C_l(u_{l_1}^pu_{l_2}^t|u_{l_1^{}}^p^{}u_{l_2^{}}^t^{}).`$ (277)
Notice that all contributions to the Fock exchange, Eqs. (272)-(274) and (277) are proportional to $`1/N`$ as one could expect from a charge distributed over $`N`$ cells. From Eq. (272) we can calculate the Coulomb self-energy associated with a conduction state $`|\stackrel{}{k},\alpha `$. Assuming $`\stackrel{}{p}=\stackrel{}{k}`$, $`\alpha =\beta `$ and taking into account the expansion (B8) and (B9), we arrive at

$`E_{si}(\stackrel{}{k},\alpha )`$ $`=`$ $`V_h+{\displaystyle \underset{\stackrel{}{K},\stackrel{}{P}}{}}{\displaystyle \underset{s_z,s_z^{}}{}}\gamma _\stackrel{}{P}^{}(\stackrel{}{k},\alpha ,s_z^{})\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z)`$ (278)
$`\times `$ $`\stackrel{}{k},\stackrel{}{P},s_z^{}|^{exc}(\stackrel{}{k},\alpha )|\stackrel{}{k},\stackrel{}{K},s_z{\displaystyle \frac{1}{N}},`$ (279)
where $`\stackrel{}{k},\stackrel{}{P},s_z^{}|^{exc}(\stackrel{}{k},\alpha )|\stackrel{}{k},\stackrel{}{K},s_z`$ is given by (272). Here we have introduced $`V_h`$ which is an electrostatic energy, associated with the electron charge $`e`$ distributed homogeneously inside the crystal. (Such homogeneous charge distribution is absent for exchange if two conduction states are different, Appendix F. It is compensated by the positive contribution of nuclei for the direct Coulomb interaction.) $`V_h`$ depends on the shape of a crystal, but decreases as $`1/N`$. As a result we observe that $`E_{si}0`$ in the limit $`N\mathrm{}`$ for any extended state.

By summing the exchange over all occupied conduction states $`|\stackrel{}{p},\beta |\stackrel{}{k},\alpha `$ we obtain a finite value of exchange for each conduction electron $`|\stackrel{}{k},\alpha `$,

$`\stackrel{}{k},\stackrel{}{P},s_z^{}|^{exc}(val)|\stackrel{}{k},\stackrel{}{K},s_z={\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{p},\beta }{}}(B_{out}^{exc}(\stackrel{}{p},\beta )`$ (280)
$`+B_I^{exc}(\stackrel{}{p},\beta )+B_S^{exc}(\stackrel{}{p},\beta )).`$ (281)
In the absence of the spin-orbit coupling, the matrix elements (272) and (281) are diagonal in spin components $`s_z=\pm 1/2`$. An important practical complication here arises due to the fact that the structural constants $`V_{e\mathrm{\hspace{0.17em}0}}`$, $`V_{e\mathrm{\Lambda }}`$ ($`V_{e\mathrm{\hspace{0.17em}0}}^{out}`$, $`V_{e\mathrm{\Lambda }}^{out}`$) depend on $`\stackrel{}{p}`$ ($`\stackrel{}{q}`$) and the calculation of them has to be repeated for each vector $`\stackrel{}{p}`$.

## VI Conclusions

We have presented a new Hartree-Fock-LAPW method for electron band structure calculations. The method combines the restricted Hartree-Fock-Roothaan approach with the crystalline basis functions in the form of linear augmented plane waves. The strategy of the full potential LAPW treatment  is adopted for calculations of the matrix elements of the direct Coulomb interactions and exchange. This is pivotal for collecting all exchange terms together including the long-range and multipole contributions.

In the framework of the FP-LAPW treatment an original technique for the solution of periodic Poisson’s equation is formulated, Sec. II. The technique takes into account the partitioning of space into two regions, inside the spheres and in the interstices, Fig. 1. In the interstitial region we expand electron densities and the potential in Fourier series and express “exchange” densities in terms of plane waves. Inside the spheres we expand densities and potentials in multipole series. Finally, we use these expansions to calculate the matrix elements of the direct Coulomb interaction (Sec. IV) and the exchange (Sec. V). The crystal field effects are considered for core electron shells and for conduction electrons. These effects are associated with the nonspherical density components of $`f`$, $`d`$ and, for noncubic symmetries, of $`p`$ electrons.  There, the crystal site symmetry is taken into account and the basis functions are adapted for the spin-orbit interaction.

The technique for solving Poisson’s equation has been applied to the face centered cubic lattice, Sec. III. We have calculated structural constants which are used to restore cubic Coulomb potentials inside a sphere from its monopole ($`l=0`$) and multipole ($`l=4,6`$) moments. We have compared our technique with the pseudo-charge-density method of Weinert, Ref.  and the two-center expansion of the Coulomb interaction, Ref. .

At present we are working on programming the formulas derived in this article. However, it is already clear that the task consists of two independent parts. First of all, one should calculate the multipole matrix elements $`c_\mathrm{\Lambda }(\lambda _1,\lambda _2)`$ of electron transitions, Eq. (166), and the other coefficients related to them (such as $`c_\mathrm{\Lambda }(\tau ,\tau ^{})`$, Eq. (159), and $`c_\mathrm{\Lambda }(\tau ,s_z;\lambda )`$, Eq. (205)). The integrations there involve only angular (and spin) parts of electronic wave functions and thus the coefficients can be tabulated and stored before the HFR self-consistent procedure. Also to this part one should add calculations of the relevant structural constants, such as $`V_{40}`$, $`V_{44}^{out}`$, $`V_{46}`$, $`V_{60}`$, and $`V_{66}^{out}`$ computed in Sec. III for the face centered cubic structure. (The constants are needed to restore the full potential for a given set of multipole moments.) These calculations depend on the type of crystal symmetry but are separated from the problems of HFR method. The second task is to program the matrix elements and all relevant procedures of the restricted HFR method. For those purposes one can start with an existing code of LAPW method and develop it on the basis of the considerations presented in this article.

###### Acknowledgements.

We thank professor K.H. Michel for numerous fruitful discussions and S. Balaban, O. Kepp and D. Kirin for helpful remarks. This work has been financially supported by the Fonds voor Wetenschappelijk Onderzoek, Vlaanderen, and by the Russian Foundation for Basic Reasearch, project No. 00-03-32968.

## A

Throughout this paper instead of complex (surface) spherical harmonics $`Y_l^m`$ we use real symmetry adapted functions  (SAFs) $`S_\mathrm{\Lambda }`$, which transform according to irreducible representations $`\mathrm{\Lambda }(l)`$ of a site symmetry group $`G`$. The composite index $`\mathrm{\Lambda }(l)`$ stands for $`(l;\mathrm{\Gamma }(l),\mu (l),k)`$ where $`\mathrm{\Gamma }(l)`$ labels the irreducible representations within the $`l`$ manifold, $`\mu (l)`$ numbers the representations that occur more than once and $`k`$ denotes the rows of a given representation. In the following we omit index $`l`$ in $`\mathrm{\Gamma }`$ and $`\mathrm{\Lambda }`$. SAFs $`S_\mathrm{\Lambda }`$ are linear combinations of $`Y_l^m`$ with the same $`l`$, i.e,

$`S_\mathrm{\Lambda }={\displaystyle \underset{m=0}{\overset{l}{}}}c_{m,c}^\mathrm{\Lambda }Y_l^{m,c}+{\displaystyle \underset{m=1}{\overset{l}{}}}c_{m,s}^\mathrm{\Lambda }Y_l^{m,s},`$ (A1)
where the coefficients $`c_{m,c}^\mathrm{\Lambda }`$ and $`c_{m,s}^\mathrm{\Lambda }`$ depend on the group $`G`$ under consideration, and $`Y_l^{m=0,c}Y_l^0`$. The coefficients $`c_{m,c,s}^\mathrm{\Lambda }`$ for different groups are quoted in Tables 2.4-2.6 of Ref. . There are $`(2l+1)`$ independent SAFs $`S_\mathrm{\Lambda }`$ belonging to the $`l`$ manifold. The real spherical harmonics are

$`Y_l^{m,c}={\displaystyle \frac{1}{\sqrt{2}}}(Y_l^m+Y_l^m)`$ (A3)
$`Y_l^{m,s}={\displaystyle \frac{1}{i\sqrt{2}}}(Y_l^mY_l^m)`$ (A4)
where $`Y_l^m`$ are taken with the phase definition of Ref. . (It is different from the definition used by Condon and Shortly. ) For a given $`l`$ in a $`(2l+1)`$ dimensional space we consider row vectors $`(Y_l^{0,c},Y_l^{1,c},\mathrm{},Y_l^{l,c},Y_l^{1,s},\mathrm{},Y_l^{l,s})Y`$ (we exclude $`Y_l^{0,s}0`$), and $`(S_\mathrm{\Lambda })S`$, and the matrix $`[{}_{c_{m,s}^\mathrm{\Lambda }}{}^{c_{m,c}^\mathrm{\Lambda }}]c`$. Then the SAFs and the spherical harmonics for a given $`l`$ are connected through an orthogonal transformation,

$`S=Yc.`$ (A5)
One can easily find the inverse transformation since

$`Y=Sc^T,`$ (A6)
where $`T`$ stands for the transpose, since $`c^1=c^T`$. It is more convenient to use $`S_\mathrm{\Lambda }`$ than $`Y_l^m`$ due to their known symmetry properties.

In the density expansion, Eq. (4), only the SAFs of $`A_{1g}`$ symmetry survive because density stays invariant under all symmetry operations of $`G`$. However, for calculations of exchange (Sec. V) there is no such simplification and the full basis set (including SAFs belonging to the other irreducible representations) should be taken into account.

We use SAFs to describe both electronic densities and wave functions. For electronic states we adopt a notation with small letters, i.e $`\lambda =(l;\mathrm{\Gamma },\nu ,k)`$. For localized electrons for conciseness we incorporate also the principal quantum number $`n`$ and write $`\tau =(n,l;\mathrm{\Gamma },\nu ,k)`$. In the latter case $`\mathrm{\Gamma }`$ refers to a double valued irreducible representations of $`G`$.

## B

Here we introduce some definitions and notations of the linear augmented plane wave method (LAPW).  The coordinate basis functions are plane waves in the interstitial region,

$`\chi _{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{R})\stackrel{}{R}|\stackrel{}{k},\stackrel{}{K}={\displaystyle \frac{1}{\sqrt{Nv}}}e^{i(\stackrel{}{k}+\stackrel{}{K})\stackrel{}{R}},`$ (B1)
and a linear combination of local atomic functions inside the spheres,

$`\stackrel{}{R}S|\stackrel{}{k},\stackrel{}{K}={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{\stackrel{}{n}}{}}e^{i\stackrel{}{k}\stackrel{}{X}(\stackrel{}{n})}{\displaystyle \underset{\lambda }{}}\varphi _\lambda ^{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{R}\stackrel{}{X}(\stackrel{}{n})),`$ (B2)
(B3)
where

$`\varphi _\lambda ^{\stackrel{}{k},\stackrel{}{K}}(\stackrel{}{r})=\left(A_\lambda (\stackrel{}{k},\stackrel{}{K})u_l(r)+B_\lambda (\stackrel{}{k},\stackrel{}{K})\dot{u}_l(r)\right)S_\lambda (\widehat{r}).`$ (B4)
(B5)
(If one uses spherical harmonics $`Y_l^m`$ instead of SAFs $`S_\lambda `$, then $`\lambda =(l,m)`$.) In order to condense notations we introduce two components ($`p`$=1,2) of the radial function $`u_l^p(r)`$,

$`u_l^1(r)u_l(r),u_l^2(r)\dot{u}_l(r)={\displaystyle \frac{u_l(r,E)}{E}}.`$ (B6)
and the corresponding to them two components $`A_\lambda ^p(\stackrel{}{k},\stackrel{}{K})`$, which are

$`A_\lambda ^1(\stackrel{}{k},\stackrel{}{K})A_\lambda (\stackrel{}{k},\stackrel{}{K}),A_\lambda ^2(\stackrel{}{k},\stackrel{}{K})B_\lambda (\stackrel{}{k},\stackrel{}{K}).`$ (B7)
The coefficients $`A_\lambda ^p`$ are obtained by requiring that the basis functions and their derivatives are continuous on the sphere boundary.

In the absence of a static magnetic field and the spin-orbit coupling, the conduction electronic states with spin projections $`s_z=\pm 1/2`$ are degenerate. (The spin-orbit coupling can be included later in the second variational treatment as described in Ref. .) The wave function of a conduction electron with the wave vector $`\stackrel{}{k}`$ and the band index $`\alpha `$ then reads

$`\stackrel{}{R}|\stackrel{}{k},\alpha ={\displaystyle \underset{s_z=1/2}{\overset{1/2}{}}}\stackrel{}{R},s_z|\stackrel{}{k},\alpha \zeta (s_z),`$ (B8)
where $`\zeta (s_z=\pm 1/2)`$ are the two spinors and

$`\stackrel{}{R},s_z|\stackrel{}{k},\alpha ={\displaystyle \underset{\stackrel{}{K}}{}}\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z)\stackrel{}{R}|\stackrel{}{k},\stackrel{}{K}.`$ (B9)
The coefficients $`\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z)`$ are found by the HF variational procedure.  Inside a sphere $`S(\stackrel{}{n})`$ the wave function is given by

$`\stackrel{}{r},s_z|\stackrel{}{k},\alpha ={\displaystyle \frac{e^{i\stackrel{}{k}\stackrel{}{X}(\stackrel{}{n})}}{\sqrt{N}}}{\displaystyle \underset{\lambda }{}}[\gamma A]_\lambda ^p(\stackrel{}{k},\alpha ,s_z)u_l^p(r)S_\lambda (\widehat{r}),`$ (B10)
(B11)
where $`\stackrel{}{r}=\stackrel{}{r}(\stackrel{}{n})`$ and summation over $`p=1,2`$ is implied. Here we have introduced the notation

$`[\gamma A]_\lambda ^p(\stackrel{}{k},\alpha ,s_z)={\displaystyle \underset{\stackrel{}{K}}{}}\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z)A_\lambda ^p(\stackrel{}{k},\stackrel{}{K}).`$ (B12)
## C

Here we derive analytical orientational wave vectors for the cubic site symmetry $`O_h`$ by employing the eigenvectors tabulated in Ref. .

The core states of $`s_{1/2}`$, $`p_{1/2}`$, $`p_{3/2}`$ and $`d_{3/2}`$ electrons remain degenerate while $`d_{5/2}`$, $`f_{5/2}`$ and $`f_{7/2}`$ are split in the cubic environment. If $`D_{5/2}`$ and $`D_{7/2}`$ are doubled valued representations of $`SO(3)`$, then the symmetry lowering is

$`D_{5/2}\mathrm{\Gamma }_7+\mathrm{\Gamma }_8,`$ (C2)
$`D_{7/2}\mathrm{\Gamma }_7+\mathrm{\Gamma }_8+\mathrm{\Gamma }_6.`$ (C3)
For two components of the doublet $`\mathrm{\Gamma }_7`$ of $`d_{5/2}`$ we have found

$`\widehat{r}|d_{5/2},(\mathrm{\Gamma }_7,1)`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{3}}}S_{2,(T2g,3)}\zeta _+`$ (C6)
$`+{\displaystyle \frac{1}{\sqrt{3}}}(iS_{2,(T2g,1)}S_{2,(T2g,2)})\zeta _{},`$
$`\widehat{r}|d_{5/2},(\mathrm{\Gamma }_7,2)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}(iS_{2,(T2g,1)}+S_{2,(T2g,2)})\zeta _+`$ (C8)
$`{\displaystyle \frac{i}{\sqrt{3}}}S_{2,(T2g,3)}\zeta _{},`$
where $`S_{2,(T2g,k)}`$ refer to the three components ($`k=13`$) of $`T_{2g}`$ symmetry ($`l=2`$) given in Table 2.6 of Ref. . The second component, Eq. (C8), is connected with the first, Eq. (C6), through the time reversal symmetry. One can easily check it by noting that the SAFs are real and by applying the following rules for the time reversal symmetry

$`ii,\zeta _+\zeta _{},\zeta _{}\zeta _+.`$ (C9)
For two components of $`\mathrm{\Gamma }_8`$ we have

$`\widehat{r}|d_{5/2},(\mathrm{\Gamma }_8,1)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{15}}}(i2S_{2,(T2g,3)}+3S_{2,(Eg,2)})\zeta _+`$ (C11)
$`+`$ $`{\displaystyle \frac{1}{\sqrt{15}}}(iS_{2,(T2g,1)}+S_{2,(T2g,2)})\zeta _{},`$ (C12)
$`\widehat{r}|d_{5/2},(\mathrm{\Gamma }_8,2)`$ $`=`$ $`\sqrt{{\displaystyle \frac{3}{5}}}S_{2,(Eg,1)}\zeta _+`$ (C13)
$``$ $`{\displaystyle \frac{1}{\sqrt{5}}}(iS_{2,(T2g,1)}+S_{2,(T2g,2)})\zeta _{}.`$ (C14)
The other two components are obtained from (C12,b) by employing the rules (C9).

We consider analogously the splitting of the $`f_{5/2}`$ ($`l=3`$) states in the cubic symmetry and obtain

$`\widehat{r}|f_{5/2},(\mathrm{\Gamma }_7,1)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{7}}}(\sqrt{3}iS_{3,A2u}+{\displaystyle \frac{2}{\sqrt{3}}}S_{3,(T2u,3)})\zeta _+`$ (C15)
$`+`$ $`{\displaystyle \frac{2}{\sqrt{21}}}(S_{3,(T2u,1)}+iS_{3,(T2u,2)})\zeta _{},`$ (C16)
and

$`\widehat{r}|f_{5/2},(\mathrm{\Gamma }_8,1)=\sqrt{{\displaystyle \frac{5}{21}}}S_{3,(T2u,3)}\zeta _++[{\displaystyle \frac{3}{2\sqrt{7}}}(S_{3,(T1u,1)}`$ (C18)
$`iS_{3,(T1u,2)})+{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{5}{21}}}(S_{3,(T2u,1)}+iS_{3,(T2u,2)})]\zeta _{},`$ (C19)
$`\widehat{r}|f_{5/2},(\mathrm{\Gamma }_8,2)=\sqrt{{\displaystyle \frac{3}{7}}}S_{3,(T1u,3)}\zeta _+{\displaystyle \frac{1}{2\sqrt{7}}}[\sqrt{3}(S_{3,(T1u,1)}`$ (C20)
$`+iS_{3,(T1u,2)})\sqrt{5}(S_{3,(T2u,1)}iS_{3,(T2u,2)})]\zeta _{}.`$ (C21)
Real SAFs $`S_{3,A2u}`$, $`S_{3,(T1u,k)}`$ and $`S_{3,(T2u,k)}`$ ($`k=13`$) are given in Table 2.6 of Ref. . The other components are obtained through the time reversal symmetry, Eq. (C9).

Finally, for $`f_{7/2}`$ core states we get

$`\widehat{r}|f_{7/2},(\mathrm{\Gamma }_6,1)={\displaystyle \frac{1}{\sqrt{3}}}(S_{3,(T1u,1)}iS_{3,(T1u,2)})\zeta _+`$ (C23)
$`{\displaystyle \frac{1}{\sqrt{3}}}S_{3,(T1u,3)}\zeta _{},`$ (C24)
$`\widehat{r}|f_{7/2},(\mathrm{\Gamma }_7,1)={\displaystyle \frac{1}{\sqrt{7}}}(2iS_{3,A2u}+S_{3,(T2u,3)})\zeta _+`$ (C25)
$`+{\displaystyle \frac{1}{\sqrt{7}}}(S_{3,(T2u,1)}+iS_{3,(T2u,2)})\zeta _{}.`$ (C26)
These are first components of $`\mathrm{\Gamma }_6`$ and $`\mathrm{\Gamma }_7`$, respectively. For the quartet $`\mathrm{\Gamma }_8`$ of $`f_{7/2}`$ we have found

$`\widehat{r}|f_{7/2},(\mathrm{\Gamma }_8,1)={\displaystyle \frac{1}{2}}[\sqrt{{\displaystyle \frac{5}{21}}}(S_{3,(T1u,1)}iS_{3,(T1u,2)})`$ (C28)
$`+{\displaystyle \frac{3}{\sqrt{7}}}(S_{3,(T2u,1)}+iS_{3,(T2u,2)})]\zeta _++\sqrt{{\displaystyle \frac{5}{21}}}S_{3,(T1u,3)}\zeta _{},`$ (C29)
$`\widehat{r}|f_{7/2},(\mathrm{\Gamma }_8,2)=\sqrt{{\displaystyle \frac{3}{7}}}S_{3,(T2u,3)}\zeta _++{\displaystyle \frac{1}{2}}[\sqrt{{\displaystyle \frac{5}{7}}}(S_{3,(T1u,1)}`$ (C30)
$`iS_{3,(T1u,2)})\sqrt{{\displaystyle \frac{3}{7}}}(S_{3,(T2u,1)}+iS_{3,(T2u,2)})]\zeta _{}.`$ (C31)
Again, the other components of $`\mathrm{\Gamma }_6`$, $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_8`$, are found by applying the time reversal transformation (C9) to (C24,b) and (C29,b).

## D

Here we describe a simple method to generate radial basis functions. We can use the spherically symmetric component of the total electrostatic potential $`V_0(\stackrel{}{r})`$ (inside a sphere $`S(\stackrel{}{n})`$) to find the radial solution $`u_{n,l}`$ and the corresponding energy $`E_{n,l}`$,

$`h_lu_{n,l}=E_{n,l}u_{n,l}.`$ (D1)
Here $`n`$ is the principal quantum number, and $`h_l`$ is the radial operator of the Schrödinger equation (Eq. (1b) of Ref. ). We assume that the solutions are confined inside the sphere $`S(\stackrel{}{n})`$ and on the sphere boundary $`r=R`$ we have

$`u_l(R)=0,{\displaystyle \frac{u_l(r)}{r}}|_R=0.`$ (D2)
The boundary conditions are complementary to (D1). Starting with (D1) and (D2) one obtains $`u_{n,l}`$ and $`E_{n,l}`$.

For a given $`l`$ we consider the radial functions $`u_{n,l}`$ which differ from each other by the principal quantum number $`n=l+1,l+2,\mathrm{}`$ . These functions correspond to the same angular dependence (specified by $`l,m_l`$) and form a complete orthonormalized set. Therefore, they can serve as a basis for any radial function $`_\tau `$ satisfying Eq. (D2). A nonrelativistic function of a localized electron is characterized by the combined index $`\tau (l,m_l)`$ and we write $`u_\eta ^\tau =u_{n,l}`$, $`i.e.`$ $`\eta =n`$. For the relativistic case one has to distinguish two solutions with the total momentum $`j=l1/2`$ and $`j=l+1/2`$. Then $`\tau =(l,j,m_j)`$ and the basis functions are $`u_\eta ^\tau `$, where again $`\eta =n`$.

For the matrix elements of the spherically symmetric potential $`V_0^\tau `$, Eq. (173), we find

$`\tau ,\eta |_0^{Coul}|\tau ,\eta ^{}`$ $`=`$ $`E_{n,J}\delta _{\eta \eta ^{}}{\displaystyle \frac{1}{4\pi }}C_0(u_\eta ^\tau u_\eta ^{}^\tau |_\tau _\tau ).`$ (D3)
Here $`E_{n,J}`$ is the eigenvalue corresponding to $`u_\eta ^\tau `$ ($`J=(l,j)`$) and the integral $`C_0(\mathrm{})`$ is given by Eq. (E2).

## E

The two-fold integral $`C_l(f_1,f_2|g_1g_2)`$ of four radial functions $`f_1(r)`$, $`f_2(r)`$, $`g_1(r)`$ and $`g_2(r)`$ is defined as

$`C_l(f_1,f_2|g_1g_2)`$ $`=`$ $`{\displaystyle _0^R}𝑑rr^2{\displaystyle _0^R}𝑑r^{}r_{}^{}{}_{}{}^{2}`$ (E2)
$`\times f_1(r)f_2(r)v_l(r,r^{})g_1(r^{})g_2(r^{})`$
where the one center multipole function is given by $`v_l(r,r^{})=\frac{4\pi }{2l+1}\frac{r_<^l}{r_>^{l+1}}`$, Eq. (16), and $`r_<`$ $`(r_>)`$ is the smaller (larger) of $`r`$ and $`r^{}`$.

## F

We consider the exchange between two extended states $`|d=|\stackrel{}{k},\alpha `$ and $`|c=|\stackrel{}{p},\beta `$. We expand both states in the LAPW basis functions, Eq. (B9). Proceeding as in section V.C, we introduce an effective “exchange” density $`\rho _{cd}(\stackrel{}{R})`$, Eq. (199), and the corresponding “exchange charge”:

$`Q_0^{cd}(R)`$ $`=`$ $`{\displaystyle \underset{l}{}}{\displaystyle \underset{\lambda (l)}{}}{\displaystyle \underset{s_z}{}}[\gamma A]_\lambda ^p(\stackrel{}{p},\beta ,s_z)[\gamma A]_\lambda ^p(\stackrel{}{k},\alpha ,s_z)N_l^p`$ (F2)
$`+𝒪(\stackrel{}{p},\beta ;\stackrel{}{k},\alpha ),`$
where

$`𝒪(\stackrel{}{p},\beta ;\stackrel{}{k},\alpha )={\displaystyle \frac{4\pi R^2}{v}}{\displaystyle \underset{\stackrel{}{K}^{}}{}}{\displaystyle \frac{j_1(|\stackrel{}{K}^{}+\stackrel{}{q}|R)}{|\stackrel{}{K}^{}+\stackrel{}{q}|}}\rho _\stackrel{}{K}^{}(\stackrel{}{p},\beta ;\stackrel{}{k},\alpha )`$ (F3)
(F4)
and

$`\rho _\stackrel{}{K}^{}(\stackrel{}{p},\beta ;\stackrel{}{k},\alpha )={\displaystyle \underset{\stackrel{}{K}}{}}{\displaystyle \underset{s_z}{}}\gamma _{\stackrel{}{K}^{}+\stackrel{}{K}}^{}(\stackrel{}{p},\beta ,s_z)\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z).`$ (F5)
On the other hand, the orthogonality relation for the two extended states is

$`\stackrel{}{p},\beta |\stackrel{}{k},\alpha =\delta _{\alpha \beta }\delta _{\stackrel{}{p}\stackrel{}{k}}`$ (F6)
$`={\displaystyle \underset{l}{}}{\displaystyle \underset{\lambda (l)}{}}{\displaystyle \underset{s_z}{}}[\gamma A]_\lambda ^p(\stackrel{}{p},\beta ,s_z)[\gamma A]_\lambda ^p(\stackrel{}{k},\alpha ,s_z)N_l^p`$ (F7)
$`+𝒪(\stackrel{}{p},\beta ;\stackrel{}{k},\alpha )+\delta _{\stackrel{}{p}\stackrel{}{k}}{\displaystyle \underset{\stackrel{}{K}}{}}{\displaystyle \underset{s_z}{}}\gamma _\stackrel{}{K}^{}(\stackrel{}{p},\beta ,s_z)\gamma _\stackrel{}{K}(\stackrel{}{k},\alpha ,s_z).`$ (F8)
Comparing it with (F2) we observe that if $`\stackrel{}{p}\stackrel{}{q}`$ then $`Q_0^{cd}(R)=0`$. However, it will not be erroneous to use this charges as they appear in HFR procedure, Eq. (245). Since Poisson’s equation is linear, their contributions will cancel in the final results. Generalizing the orthonormality relation (F8) for the case of few atoms one can show that it leads to

$`{\displaystyle \underset{i}{}}Q_i^{cd}=0,\stackrel{}{p}\stackrel{}{k}.`$ (F9)
We conclude that the “exchange” charges $`Q_i^{cd}`$ are not necessarily zero. (Here $`i`$ labels different atoms in the unit cell.) The orthogonality relation (F9) ensures that there is no Coulomb divergence associated with the uniform component of the “exchange density”.