Spin-polarized relativistic calculations in linear augmented-Slater-type-orbital method

Let us establish a fundamental one-electronic equation, before the DLASTO method will be formulated. Here it is assumed that the spin-polarized relativistic density-functional is taken up to describe the electronic states in magnetic materials or atoms we are interested in. For such the systems the total energy in an external field ($V\left(\boldsymbol{r}\right)$, A ( r )) is defined as a functional of four-dimensional current (n( r ), J ( r )), and that the ground state is derived by taking a minimum of the functional so as to obtain a relativistic one-electronic equation [36]. In order to derive a numerical expression, the Gordon decomposition of the current is applied to this equation [37], thereby giving a useful equation of the 'spin-only' type [38–40] in units of Ry as,

$$\begin{equation}\left[c\boldsymbol{\alpha }\cdot \boldsymbol{p}+(\enspace \beta -I)\frac{{c}^{2}}{2}+V(\boldsymbol{r})+\beta \mathbf{\Sigma }\cdot \boldsymbol{B}(\boldsymbol{r})\right]{{\Phi}}_{\lambda }(\boldsymbol{k},\boldsymbol{r})={E}_{\lambda }(\boldsymbol{k}){{\Phi}}_{\lambda }(\boldsymbol{k},\boldsymbol{r}),\end{equation} \tag{ 1 }$$

where c is the light velocity and p means the momentum operator (≡−i∇). This three operators, α , β and I represent the standard Dirac 4 × 4 matrices [27], and Σ signifies the 4 × 4 Pauli matrices. Then E_λ ( k ) is a band energy identified with band index λ and k vector in the first Brillouin zone, and Φ_λ ( k , r ) denotes the corresponding Bloch function.

In this formulation, a crystal structure is assumed to consist of n_a different groups of atoms as a crystalline basis in the unit cell, and the basis contains n_p atoms in the pth atomic group. The νth atom of the pth group (the pνth atom) is directed by vector ${\boldsymbol{\tau }}_{\nu }^{p}$ from an origin in the unit cell. As in the usual band theories, a real space of the unit cell is divided into the sphere region inside the spheres around each of the atoms and the interstitial region outside the spheres. Here a radius of the sphere denotes a_p for the atoms of the pth group.

In (1), V( r ) and B ( r ) show the scalar potential and the magnetic field, respectively, each of which is defined by,

$$\begin{equation}V(\boldsymbol{r})=-\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{\nu =1}^{{n}_{p}}\sum\limits _{\boldsymbol{R}}\frac{2{Z}_{p}}{\vert \boldsymbol{r}-{\boldsymbol{\tau }}_{\nu }^{p}+\boldsymbol{R}\vert }+\int \frac{2n({\boldsymbol{r}}^{\prime })}{\vert \boldsymbol{r}-{\boldsymbol{r}}^{\prime }\vert }\enspace \mathrm{d}{\boldsymbol{r}}^{\prime }+\frac{\partial \left[n(\boldsymbol{r}){\varepsilon }_{\mathrm{X}\mathrm{C}}(\boldsymbol{r})\right]}{\partial n(\boldsymbol{r})},\end{equation} \tag{ 2 }$$

using the translation vector R and

$$\begin{equation}\mathbf{B}(\boldsymbol{r})=\frac{\partial \left[n(\boldsymbol{r}){\varepsilon }_{\mathrm{X}\mathrm{C}}(\boldsymbol{r})\right]}{\partial \mathbf{m}(\boldsymbol{r})},\end{equation} \tag{ 3 }$$

with the charge density

$$\begin{equation}n(\boldsymbol{r})=\sum\limits _{\boldsymbol{k}\lambda }{{\Phi}}_{\lambda }^{{\dagger}}(\boldsymbol{k},\boldsymbol{r}){{\Phi}}_{\lambda }(\boldsymbol{k},\boldsymbol{r})\theta ({E}_{\mathrm{F}}-{E}_{\lambda }(\boldsymbol{k})),\end{equation} \tag{ 4 }$$

and the magnetization

$$\begin{equation}\mathbf{m}(\boldsymbol{r})=\sum\limits _{\boldsymbol{k}\lambda }{{\Phi}}_{\lambda }^{{\dagger}}(\boldsymbol{k},\boldsymbol{r})\beta \mathbf{\Sigma }{{\Phi}}_{\lambda }(\boldsymbol{k},\boldsymbol{r})\theta ({E}_{\mathrm{F}}-{E}_{\lambda }(\boldsymbol{k})),\end{equation} \tag{ 5 }$$

where θ(x) is the step function and ɛ_XC( r ) represents the exchange–correlation (XC) energy per a particle. In (4) and (5), E_F stands for the Fermi energy, so that the sum or the integration is done over all the occupied states of λ and k in the first Brillouin zone. The first term described in the right side of (2) is the Coulomb potential between electrons and nuclei with atomic number Z_p in the pth atomic group and the second term denotes the Coulomb potential between electrons. The third term of (2) and the expression of (3) indicate the scalar and the spin-dependent part of the XC potentials, respectively and the latter represents an internal magnetic field.

As in the full potential method, the scalar potential V(r) and the magnetic field B( r ) are expressed by a Fourier expansion in the interstitial region, while in the sphere region they are expanded by a spherical-harmonics. By taking account of the crystal symmetry, the interstitial scalar potential and magnetic field are explicitly described as

$$\begin{equation}V(\boldsymbol{r})=\sum\limits _{\boldsymbol{G}}\;~{v}(\boldsymbol{G})\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\boldsymbol{G}\cdot \boldsymbol{r})=\sum\limits _{{G}_{s}}~{V}({G}_{s}){\varphi }_{s}(\boldsymbol{r}),\end{equation} \tag{ 6 }$$

and

$$\begin{equation}\mathbf{B}(\boldsymbol{r})=\sum\limits _{{G}_{s}}~{\mathbf{B}}({G}_{s}){\varphi }_{s}(\boldsymbol{r}),\end{equation} \tag{ 7 }$$

where G denotes a reciprocal-lattice vector, φ_s (r) are symmetrized plane-wave functions, or stars [14], which are invariant under rotational operators in the crystalline space group and G_s means a reciprocal modulus of generators to create the sth star φ_s (r). In the sphere region we can write down the scalar potential by

$$\begin{equation}V(\boldsymbol{r})=\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{t}{V}_{t}^{p}(r){K}_{t}^{p}(\hat{\boldsymbol{r}}),\end{equation} \tag{ 8 }$$

and the magnetic field as

$$\begin{equation}\mathbf{B}(\boldsymbol{r})=\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{t}{B}_{t}^{p}(r){K}_{t}^{p}(\hat{\boldsymbol{r}}),\end{equation} \tag{ 9 }$$

where ${K}_{t}^{p}(\hat{\boldsymbol{r}})$ means normalized lattice harmonics with the site symmetry of pth atomic group [14] that keep all the atomic positions invariant. The hat-notation of $\hat{\boldsymbol{r}}$ stands for an angular part in a polar coordinate and it is henceforth used in formulation. The explicit definition of ${K}_{t}^{p}(\hat{\boldsymbol{r}})$ used here is given in (A.11) of appendix A.

The Bloch function Φ( k , r ) is expanded by DLASTO basis function ${{\Psi}}_{n\ell s\mu }^{p\nu }(\boldsymbol{k},\boldsymbol{r})$ as follows,

$$\begin{equation}{\Phi}(\boldsymbol{k},\boldsymbol{r})=\sum\limits _{n\ell s\mu }\sum\limits _{p\nu }{U}_{n\ell s\mu }^{p\nu }{{\Psi}}_{n\ell s\mu }^{p\nu }(\boldsymbol{k},\boldsymbol{r})=\sum\limits _{N}{U}_{N}{{\Psi}}_{N}(\boldsymbol{k},\boldsymbol{r}),\end{equation} \tag{ 10 }$$

where the index N is an abbreviation for {nℓsμpν} and the band index λ is omitted here, though it is present in the expansion coefficients ${U}_{n\ell s\mu }^{p\nu }$.

In a non-relativistic case, the basis function in the interstitial region can be described using Fourier components of the Bloch sum of STO [17], that are atomic orbital functions simplified with principal quantum number n, angular momentum ℓ and screening factor ζ [10], expressed as,

$$\begin{equation}{\varphi }_{nl}(\boldsymbol{r})={r}^{n-1}\enspace \mathrm{e}\mathrm{x}\mathrm{p}(-\zeta r){Y}_{\ell }^{m}(\hat{\boldsymbol{r}}).\end{equation} \tag{ 11 }$$

As some analytically-derived results, a Dirac-type relativistic generalization to the interstitial basis function can be simply done to replace the STO or the plane waves with any four-component type of wave function. As a choice, here, we define the interstitial basis by use of the Fourier transformation of φ_nl ( r ) with relativistic plane waves [27], because it is highly compatibility to the DLAPW method [21]. If it is centered at the νth atom associated with the pth group and if its quantum state of orbital is in the principal n, the angular ℓ, the total angular J = ℓ + s and its z-component μ, then the basis function can be defined within the order of 1/c by,

$$\begin{equation}{{\Psi}}_{N}(\boldsymbol{k},\boldsymbol{r})=\frac{1}{\sqrt{{\Omega}}}\sum\limits _{{\boldsymbol{k}}_{j}}\sum\limits _{m={\pm}\frac{1}{2}}{S}_{N}({\boldsymbol{k}}_{j},m){{\Psi}}_{\text{PW}}({\boldsymbol{k}}_{j},m;\boldsymbol{r}),\end{equation} \tag{ 12 }$$

where Ψ_PW( k _j , m; r ) are the relativistic plan waves, represented explicitly in units of Ry as,

$$\begin{equation}{{\Psi}}_{\text{PW}}(\boldsymbol{k},m;\boldsymbol{r})=\sqrt{\frac{2{k}^{\mathrm{o}}+{c}^{2}}{4{k}^{\mathrm{o}}}}\left(\begin{matrix}\hfill \chi (m)\hfill \\ \hfill \frac{2c\boldsymbol{\sigma }\cdot \boldsymbol{k}}{2{k}^{\mathrm{o}}+{c}^{2}}\chi (m)\hfill \end{matrix}\right)\mathrm{e}\mathrm{x}\mathrm{p}[\mathrm{i}\boldsymbol{k}\cdot \boldsymbol{r}]\enspace .\end{equation} \tag{ 13 }$$

Then S_N ( k , m) are the envelop functions of Ψ_PW( k , m; r ),

$$\begin{equation}{S}_{N}(\boldsymbol{k},m)={S}_{n\ell s\mu }^{p\nu }(\boldsymbol{k},m)=\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}\boldsymbol{k}\cdot {\boldsymbol{\tau }}_{\nu }^{p}){~{\varphi }}_{n\ell s\mu }^{p}(\boldsymbol{k},m),\end{equation} \tag{ 14 }$$

where Ω is the volume of unit cell, and k _j is defined using the reciprocal-lattice vector G _j as k _j = k + G _j . In (13), k^o is the relativistic energy

$$\begin{equation}2{k}^{\mathrm{o}}=c\sqrt{{(2k)}^{2}+{c}^{2}},\end{equation} \tag{ 15 }$$

with a wavevector modulus k. In (13), σ are the usual Pauli matrices and χ(m) signifies the spinor function for $m={\pm}\frac{1}{2}$. In (14) here we use an index $s={\pm}\frac{1}{2}$ instead of J, and we should note that, for ℓ = 0, only $s=\frac{1}{2}$ is available in the $s={\pm}\frac{1}{2}$ formulation.

In (14), ${~{\varphi }}_{n\ell s\mu }^{p}(\boldsymbol{k},m)$ corresponds to the main term for the Fourier components of relativistic STO, which is described as,

$$\begin{equation}{~{\varphi }}_{n\ell s\mu }^{p}(\boldsymbol{k},m)={~{\phi }}_{n\ell }^{p}(k,{\zeta }_{p}){\mathcal{Y}}_{\ell s}^{\mu }(\hat{\boldsymbol{k}},m),\end{equation} \tag{ 16 }$$

with the Fourier component of radial STO for the pth group,

$$\begin{equation}{~{\phi }}_{n\ell }^{p}(k,\zeta )={\int }_{0}^{\infty }{r}^{n+1}\enspace \mathrm{e}\mathrm{x}\mathrm{p}(-\zeta r){j}_{\ell }(kr)\mathrm{d}r={(-1)}^{n}{\left[\frac{\partial }{\partial \zeta }\right]}^{n-\ell }{\left[\frac{1}{\zeta }\frac{\partial }{\partial \zeta }\right]}^{\ell }\frac{{k}^{\ell }}{{k}^{2}+{\zeta }^{2}},\end{equation} \tag{ 17 }$$

and the relativistic angular dependent term,

$$\begin{equation}{\mathcal{Y}}_{\ell s}^{\mu }(\hat{k},m)={C}_{sm}^{\ell ;\mu }{Y}_{\ell }^{\mu -m}(\hat{\boldsymbol{k}}),\end{equation} \tag{ 18 }$$

where j_ℓ (x) is the spherical Bessel function of ℓ order. Practically the numerical calculation of ${~{\phi }}_{n\ell }^{p}(k,\zeta )$ in (17) can be directly treated with terms of the Gauss hypergeometric function [41]. Then ${C}_{sm}^{\ell ;\mu }$ shows the Clebsch–Gordan (CG) coefficients, which are defined with the parameter ${u}_{\mu }^{\ell }=\mu {(\ell +1/2)}^{-1}$ as the matrices C ^ℓμ ,

$$\begin{equation}{\boldsymbol{C}}^{\ell ;\mu }=\left(\begin{matrix}\hfill {C}_{\frac{1}{2}\frac{1}{2}}^{\ell ;\mu }\hfill & \hfill {C}_{\frac{1}{2}-\frac{1}{2}}^{\ell ;\mu }\hfill \\ \hfill {C}_{-\frac{1}{2}\frac{1}{2}}^{\ell ;\mu }\hfill & \hfill {C}_{-\frac{1}{2}-\frac{1}{2}}^{\ell ;\mu }\hfill \end{matrix}\right)=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill \sqrt{1+{u}_{\mu }^{\ell }}\hfill & \hfill \sqrt{1-{u}_{\mu }^{\ell }}\hfill \\ \hfill -\sqrt{1-{u}_{\mu }^{\ell }}\hfill & \hfill \sqrt{1+{u}_{\mu }^{\ell }}\hfill \end{matrix}\right).\end{equation} \tag{ 19 }$$

The STO screening factor ζ_p in (16) is determined from the boundary condition on the sphere surrounding any atom of the pth group by a common way of the linear band theories D_ℓ = −ℓ − 1 [13] and ζ_p is derived from the following expression [17],

$$\begin{equation}{\zeta }_{p}=\frac{n+\ell }{{a}_{p}}.\end{equation} \tag{ 20 }$$

Inside the sphere at the pνth atom, the basis Ψ_N ( k , r ) in (10) is augmented by a linear combination of wave function ${\psi }_{\ell \mu }^{\alpha p}(\boldsymbol{r})$ and its corresponding energy-derivative wave function ${\dot {\psi }}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha p}(\boldsymbol{r})=\partial {\psi }_{\ell \mu }^{\alpha p}(\boldsymbol{r})/\partial \varepsilon$, and is described by,

$$\begin{equation}{{\Psi}}_{N}(\boldsymbol{k},\boldsymbol{r})=\sum\limits _{{\ell }^{\prime }{\mu }^{\prime }}\sum\limits _{\alpha ={\pm}\frac{1}{2}}\left[{B}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k}){\psi }_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha p}(\boldsymbol{r}-{\boldsymbol{\tau }}_{\nu }^{p})+{\dot {B}}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k}){\dot {\psi }}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha p}(\boldsymbol{r}-{\boldsymbol{\tau }}_{\nu }^{p})\right].\end{equation} \tag{ 21 }$$

The wave function ${\psi }_{\ell \mu }^{\alpha p}(\boldsymbol{r})$, which is a solution of (1) around the pν atom with a spherically symmetric potentials ${V}^{p}(r)=\frac{1}{\sqrt{4\pi }}{V}_{0}^{p}(r)$ of (8) and ${B}^{p}(r)=\frac{1}{\sqrt{4\pi }}{B}_{0}^{p}(r)$ of (9) along spin-quantized z-axis, is expressed as,

$$\begin{equation}{\psi }_{\ell \mu }^{\alpha p}(\boldsymbol{r})=\frac{1}{r}\sum\limits _{s={\pm}\frac{1}{2}}\left(\begin{matrix}\hfill {g}_{s}^{\ell \mu ;\alpha p}(r){\chi }_{s}^{\ell ;\mu }(\hat{\boldsymbol{r}})\hfill \\ \hfill \mathrm{i}{f}_{s}^{\ell \mu ;\alpha p}(r){\chi }_{-s}^{\ell +2s;\mu }(\hat{\boldsymbol{r}})\hfill \end{matrix}\right),\end{equation} \tag{ 22 }$$

where ${g}_{s}^{\ell \mu ;\alpha p}(r)$ and ${f}_{s}^{\ell \mu ;\alpha p}(r)$ are called the large and the small components of the relativistic radial functions, respectively. It should be noted that (22) is only validated by an approximation that breaks the orbital coupling with ℓ ± 2 [24, 32, 33]. The function ${\chi }_{s}^{\ell ;\mu }(\hat{\mathbf{r}})$ represents a normalized spin-angular function,

$$\begin{equation}{\chi }_{s}^{\ell ;\mu }(\hat{\boldsymbol{r}})=\sum\limits _{m={\pm}1/2}{C}_{sm}^{\ell ;\mu }{Y}_{l}^{\mu -m}(\hat{\boldsymbol{r}})\chi (m),\end{equation} \tag{ 23 }$$

with the CG coefficients ${C}_{sm}^{\ell ;\mu }$ in (19) and the spinor function χ(m).

In a spin-unpolarized case, α equals to $s={\pm}\frac{1}{2}$ due to diagonalizing the states with j = ℓ + α and then ${g}_{\alpha }^{\ell \mu ;\alpha p}(r)$ and ${f}_{\alpha }^{\ell \mu ;\alpha p}(r)$ become a solution of the usual Dirac equation. If the radial magnetic field B(r) is nonzero, the spin polarization has ${g}_{s}^{\ell \mu ;\alpha p}(r)$ and ${f}_{s}^{\ell \mu ;\alpha p}(r)$ coupled between the j = ℓ + s states with $s={\pm}\frac{1}{2}$. With the parameter ${u}_{\mu }^{\ell }=\mu {(\ell +1/2)}^{-1}$ used in the CG coefficients (19), the simultaneous differential equations to (22) can be written down in a compact form as follows,

$$\begin{equation}\left(\frac{\mathrm{d}}{\mathrm{d}r}+\frac{{\kappa }_{s}^{\ell }}{r}\right){g}_{s}^{\ell \mu ;\alpha p}(r)-\left[1+\frac{1}{{c}^{2}}({E}_{\ell \mu }^{\alpha p}-{V}^{p}(r)-2s{B}^{p}(r){u}_{\mu }^{\ell +2s})\right]c{f}_{s}^{\ell \mu ;\alpha p}(r)=0,\end{equation} \tag{ 24 }$$

$$\begin{align}\left(\frac{\mathrm{d}}{\mathrm{d}r}-\frac{{\kappa }_{s}^{\ell }}{r}\right)c{f}_{s}^{\ell \mu ;\alpha p}(r)+\left({E}_{\ell \mu }^{\alpha p}-{V}^{p}(r)-2s{B}^{p}(r){u}_{\mu }^{\ell }\right){g}_{s}^{\ell \mu ;\alpha p}(r)\\ +{B}^{p}(r)\sqrt{1-{({u}_{\mu }^{\ell })}^{2}}{g}_{-s}^{\ell \mu ;\alpha p}(r)=0,\end{align} \tag{ 25 }$$

using ${\kappa }_{s}^{\ell }=-2s(\ell +s+\frac{1}{2})$ for $s={\pm}\frac{1}{2}$ [24, 26]. By the presence of the last term in (25), we have to handle simultaneously the four differential equations for $s={\pm}\frac{1}{2}$ of (24) and (25). Here they are called SPCD equations and a computational algorithm to solve the SPCD equation is described in details elsewhere [26].

In a similar way the energy-derivative wave function ${\dot {\psi }}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha p}(\boldsymbol{r})$ is defined by the style of (22), in which just only ${g}_{s}^{\ell \mu ;\alpha p}(r)$ and ${f}_{s}^{\ell \mu ;\alpha p}(r)$ are replaced to ${\dot {g}}_{s}^{\ell \mu ;\alpha p}(r)$ and ${\dot {f}}_{s}^{\ell \mu ;\alpha p}(r)$, respectively, and the corresponding SPCD equations for ${\dot {g}}_{s}^{\ell \mu ;\alpha p}(r)$ and ${\dot {f}}_{s}^{\ell \mu ;\alpha p}(r)$ are given as,

$$\begin{align}\left(\frac{\mathrm{d}}{\mathrm{d}r}+\frac{{\kappa }_{s}^{\ell }}{r}\right){\dot {g}}_{s}^{\ell \mu ;\alpha p}(r)-\left[1+\frac{1}{{c}^{2}}({E}_{\ell \mu }^{\alpha p}-{V}^{p}(r)-2s{B}^{p}(r){u}_{\mu }^{\ell +2s})\right]c{\dot {f}}_{s}^{\ell \mu ;\alpha p}(r)=\frac{1}{c}{f}_{s}^{\ell \mu ;\alpha p}(r),\end{align} \tag{ 26 }$$

$$\begin{align}\left(\frac{\mathrm{d}}{\mathrm{d}r}-\frac{{\kappa }_{s}^{\ell }}{r}\right)c{\dot {f}}_{s}^{\ell \mu ;\alpha p}(r)+\left({E}_{\ell \mu }^{\alpha p}-{V}^{p}(r)-2s{B}^{p}(r){u}_{\mu }^{\ell }\right){\dot {g}}_{s}^{\ell \mu ;\alpha p}(r)\\ \enspace +B{(r)}^{p}\sqrt{1-{({u}_{\mu }^{\ell })}^{2}}{\dot {g}}_{-s}^{\ell \mu ;\alpha p}(r)=-{g}_{s}^{\ell \mu ;\alpha p}(r).\end{align} \tag{ 27 }$$

In (24)–(27), ${E}_{\ell \mu }^{\alpha p}$ is a fixed energy with respect to the (ℓνα) state in the pth atomic group [21].

The two coefficients of ${B}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})$ and ${\dot {B}}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})$ in (21) are determined from the condition that the interstitial (12) and the sphere (21) of the basis set are continuous on the spheres. For N = {nℓsμpν}, the coefficients are obtained by solving the following simultaneous equations,

$$\begin{equation}\sum\limits _{\alpha ={\pm}\frac{1}{2}}\left[{B}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k}){g}_{{s}^{\prime }}^{{\ell }^{\prime }{\mu }^{\prime };\alpha p}({a}_{p})+{\dot {B}}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k}){\dot {g}}_{{s}^{\prime }}^{{\ell }^{\prime }{\mu }^{\prime };\alpha p}({a}_{p})\right]={G}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}(N;\boldsymbol{k}),\end{equation} \tag{ 28 }$$

$$\begin{equation}\sum\limits _{\alpha ={\pm}\frac{1}{2}}\left[{B}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})c{f}_{{s}^{\prime }}^{{\ell }^{\prime }{\mu }^{\prime };\alpha {p}^{\prime }}({a}_{{p}^{\prime }})+{\dot {B}}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})c{\dot {f}}_{{s}^{\prime }}^{{\ell }^{\prime }{\mu }^{\prime };\alpha {p}^{\prime }}({a}_{{p}^{\prime }})\right]=-2{s}^{\prime }{F}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}(N;\boldsymbol{k}),\end{equation} \tag{ 29 }$$

using

$$\begin{equation}{G}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}(N;\boldsymbol{k})=\sum\limits _{{\boldsymbol{k}}_{j}}\sum\limits _{{m}^{\prime }m}{P}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}^{p{\ast}}({\boldsymbol{k}}_{j},{m}^{\prime }){D}_{{m}^{\prime }m}({\theta }_{\nu }^{p},{\phi }_{\nu }^{p}){~{\varphi }}_{n\ell s\mu }^{p}({\boldsymbol{k}}_{j},m),\end{equation} \tag{ 30 }$$

and

$$\begin{equation}{F}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}(N;\boldsymbol{k})=\sum\limits _{{\boldsymbol{k}}_{j}}\sum\limits _{{m}^{\prime }m}{Q}_{{\ell }^{\prime }{s}^{\prime }{\mu }^{\prime }}^{p{\ast}}({\boldsymbol{k}}_{j},{m}^{\prime }){D}_{{m}^{\prime }m}({\theta }_{\nu }^{p},{\phi }_{\nu }^{p}){~{\varphi }}_{n\ell s\mu }^{p}({\boldsymbol{k}}_{j},m),\end{equation} \tag{ 31 }$$

with the spherical-coordinate coefficients of

$$\begin{equation}{P}_{\ell s\mu }^{p}(\boldsymbol{k},m)=\frac{4\pi {(-\mathrm{i})}^{\ell }{a}_{p}}{\sqrt{{\Omega}}}{j}_{\ell }(k{a}_{p}){\mathcal{Y}}_{\ell \mu s}(\hat{\boldsymbol{k}},m),\end{equation} \tag{ 32 }$$

and

$$\begin{equation}{Q}_{\ell s\mu }^{p}(\boldsymbol{k},m)=\frac{4\pi {(-\mathrm{i})}^{\ell }{a}_{p}}{\sqrt{{\Omega}}}k{j}_{\ell +2s}(k{a}_{p}){\mathcal{Y}}_{\ell \mu s}(\hat{\boldsymbol{k}},m).\end{equation} \tag{ 33 }$$

In (30) and (31), ${D}_{{m}^{\prime }m}({\theta }_{\nu }^{p},{\phi }_{\nu }^{p})$ is the m'm element of spin-rotational matrix,

$$\begin{equation}\mathbf{D}(\theta ,\phi )=\left(\begin{matrix}\hfill \mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\frac{1}{2}\phi \right)\mathrm{cos}\left(\frac{1}{2}\theta \right)\hfill & \hfill \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}\frac{1}{2}\phi \right)\mathrm{sin}\left(\frac{1}{2}\theta \right)\hfill \\ \hfill -\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\frac{1}{2}\phi \right)\mathrm{sin}\left(\frac{1}{2}\theta \right)\hfill & \hfill \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}\frac{1}{2}\phi \right)\mathrm{cos}\left(\frac{1}{2}\theta \right)\hfill \end{matrix}\right),\end{equation} \tag{ 34 }$$

with respect to the Euler angles, ${\theta }_{\nu }^{p}$ and ${\phi }_{\nu }^{p}$, at the pνth atom. In non-collinear magnets that a local quantization axis of the spin in the sphere region is not along the z axis of a crystalline global coordinate, we can self-consistently calculate the local spin moment at the pνth atom by rotating the local axis by the angles ${\theta }_{\nu }^{p}$ and ${\phi }_{\nu }^{p}$ from the global coordinate. Using a relativistic density-matrix, the Euler angles, ${\theta }_{\nu }^{p}$ and ${\phi }_{\nu }^{p}$ can be determined also from a self-consistent calculation, but here we describe only the expressions to compute a fixed magnetic structure.

In the derivation of (28)–(33), we use the same approximation of k_j ≪ c as in the relativistic APW method [35]. It is to be noted that, although N contains p and ν, G_ℓ's'μ'(N; k ) in (30) and F_ℓ's'μ'(N; k ) in (31) are independent of the ${\boldsymbol{\tau }}_{\nu }^{p}$ vectors. Accordingly ${B}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})$ and ${\dot {B}}_{{\ell }^{\prime }{\mu }^{\prime }}^{\alpha }(N;\boldsymbol{k})$ have the same values for all the atoms of the pth group.

The band energies E_λ ( k ) and the corresponding eigenvectors ${U}_{N}^{\lambda }$ for the band index λ are obtained by solving the set of equation,

$$\begin{equation}\sum\limits _{{N}^{\prime }}\left[{H}_{N{N}^{\prime }}-E(\boldsymbol{k}){O}_{N{N}^{\prime }}\right]{U}_{{N}^{\prime }}=0\end{equation} \tag{ 35 }$$

with the Hamiltonian matrix elements

$$\begin{equation}{H}_{N{N}^{\prime }}={\int }_{{\Omega}}\mathrm{d}\boldsymbol{r}{{\Psi}}_{N}^{{\dagger}}(\boldsymbol{k},\boldsymbol{r})\hat{H}{{\Psi}}_{{N}^{\prime }}(\boldsymbol{k},\boldsymbol{r})={H}_{N{N}^{\prime }}^{\text{MT}}+{H}_{N{N}^{\prime }}^{\text{NMT}},\end{equation} \tag{ 36 }$$

and the overlapping matrix elements

$$\begin{equation}{O}_{N{N}^{\prime }}={\int }_{{\Omega}}\mathrm{d}\boldsymbol{r}{{\Psi}}_{N}^{{\dagger}}(\boldsymbol{k},\boldsymbol{r}){{\Psi}}_{{N}^{\prime }}(\boldsymbol{k},\boldsymbol{r}).\end{equation} \tag{ 37 }$$

The Hamiltonian matrix elements H_NN' in (36) are separated into a muffin-tin (MT) term ${H}_{N{N}^{\prime }}^{\text{MT}}$ and a non-muffin-tin (NMT) term ${H}_{N{N}^{\prime }}^{\text{NMT}}$. The explicit expressions of ${H}_{N{N}^{\prime }}^{\text{MT}}$ are casted in (A.1) of appendix A, together with O_NN' in (A.7). The NMT term ${H}_{N{N}^{\prime }}^{\text{NMT}}$ is expressed by

$$\begin{equation}{H}_{N{N}^{\prime }}^{\text{NMT}}={\int }_{{\Omega}}\mathrm{d}\boldsymbol{r}{{\Psi}}_{N}^{{\dagger}}(\boldsymbol{k},\boldsymbol{r})\left({V}^{\text{NMT}}(\boldsymbol{r})+\beta \mathbf{\Sigma }\cdot {\boldsymbol{B}}^{\text{NMT}}(\boldsymbol{r})\right){{\Psi}}_{{N}^{\prime }}(\boldsymbol{k},\boldsymbol{r}),\end{equation} \tag{ 38 }$$

where V^NMT( r ) and B ^NMT( r ) means the non-spherical terms in the scalar potential V( r ) of (2) and the magnetic field B ( r ) of (3), respectively. In fact, V^NMT( r ) is obtainable from extracting the MT terms of $~{V}(0)$ and V^p (r) from V( r ), and in the same manner B ^NMT( r ) is given by subtracting $~{B}(0)$ and B^p (r) from B ( r ). The details of ${H}_{N{N}^{\prime }}^{\text{NMT}}$ are in (A.8)–(A.18).

Once the energies and the eigenvectors for the band index λ are calculated from the general eigenvalue solutions of (35), the resulting eigenvectors ${U}_{N}^{\lambda }$ are transformed into ${~{U}}^{\lambda }({\boldsymbol{k}}_{j},m)$ with a conversion equation

$$\begin{equation}{~{U}}_{\lambda }^{\boldsymbol{k}}({\boldsymbol{G}}_{j},m)=\sum\limits _{N}{U}_{N}^{\lambda }{S}_{N}({\boldsymbol{k}}_{j},m).\end{equation} \tag{ 39 }$$

Thereby the Bloch function in the interstitial region is rewritten from (10) and (12) as

$$\begin{equation}{{\Phi}}^{\lambda }(\boldsymbol{k},\boldsymbol{r})=\frac{1}{\sqrt{{\Omega}}}\sum\limits _{j}\sum\limits _{m={\pm}\frac{1}{2}}{~{U}}_{\lambda }^{\boldsymbol{k}}({\boldsymbol{G}}_{j},m){{\Psi}}_{\text{PW}}({\boldsymbol{k}}_{j},m;\boldsymbol{r}).\end{equation} \tag{ 40 }$$

That is really the LAPW interstitial Bloch function, so that we can handle a prescription of a full-potential LAPW method [42] to develop the full potential scheme.

In obtaining the charge density (4) and the magnetization (5), the sums or the integrations over the occupied states are needed to perform, and to be more precise the tetrahedron method [43] is used on the assumption that any metallic system will be calculated here. For a semiconductor or an insulator the special point method is favorable to use [44]. From the tetrahedron method or the special point method, if the weights W^λ ( k ) are set out for the occupied states of λ and k in the irreducible wedge of the first Brillouin zone, the interstitial charge density n( r ) and magnetization m( r ) along the z axis of the global coordinate are calculated as follows,

$$\begin{equation}n(\boldsymbol{r})=\rho \left(\boldsymbol{r},+\frac{1}{2}\right)+\rho \left(\boldsymbol{r},-\frac{1}{2}\right)=\sum\limits _{{G}_{s}}~{n}({G}_{s}){\varphi }_{s}(\boldsymbol{r}),\end{equation} \tag{ 41 }$$

$$\begin{equation}m(\boldsymbol{r})=\rho \left(\boldsymbol{r},+\frac{1}{2}\right)-\rho \left(\boldsymbol{r},-\frac{1}{2}\right)=\sum\limits _{{G}_{s}}~{m}({G}_{s}){\varphi }_{s}(\boldsymbol{r}),\end{equation} \tag{ 42 }$$

using a spin-dependent density with $m={\pm}\frac{1}{2}$,

$$\begin{equation}\rho (\boldsymbol{r},m)=\frac{1}{{\Omega}}\sum\limits _{\boldsymbol{k}\lambda }{W}^{\lambda }(\boldsymbol{k})\sum\limits _{ij}{~{U}}_{\lambda }^{\boldsymbol{k}{\ast}}({\boldsymbol{G}}_{i},n){~{U}}_{\lambda }^{\boldsymbol{k}}({\boldsymbol{G}}_{j},n)\mathrm{e}\mathrm{x}\mathrm{p}[\mathrm{i}({\boldsymbol{G}}_{j}-{\boldsymbol{G}}_{i})\cdot \boldsymbol{r}]\enspace \end{equation} \tag{ 43 }$$

$$\begin{equation}=\sum\limits _{{G}_{s}}~{\rho }({G}_{s},n){\varphi }_{s}(\boldsymbol{r}).\end{equation} \tag{ 44 }$$

The right term of (43) is derived approximately from k_j ≪ c and then it is transformed to (44) by symmetrization to the stars φ_s ( r ).

In an analogous way, the charge density in the sphere region is settled from (4) and (21) with the normalized lattice harmonics ${K}_{t}^{p}(\hat{\boldsymbol{r}})$ as

$$\begin{equation}n(\boldsymbol{r})=\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{t}\;{n}_{t}^{p}(r){K}_{t}^{p}(\hat{\boldsymbol{r}}),\end{equation} \tag{ 45 }$$

where ${n}_{t}^{p}(r)$ is a radial contribution to the pth atomic group and the details are shown in (B.1)–(B.4) of appendix B. Similarly the magnetization along the z axis of the local coordinate in the sphere region is described by

$$\begin{equation}m(\boldsymbol{r})=\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{t}{m}_{t}^{p}(r){K}_{t}^{p}(\hat{\boldsymbol{r}}),\end{equation} \tag{ 46 }$$

where ${m}_{t}^{p}(r)$ is a radial contribution of magnetization in the pth atomic group, of which the expressions are written down in (B.5) and (B.6).

The interstitial and sphere scalar potentials V( r ) of (6) and (8) are constructed from the densities of (41), (42), (45) and (46). As a routine application to the full-potential scheme we need a special procedure to yield the interstitial Coulomb potentials in the first and second terms of the right side in (2), that is known as the pseudo-charge method [42, 45]. In order to achieve a well convergence within a finite number of expansion, the interstitial charge density of (41) including the nuclear charges should be replaced to a pseudo-density, which varies with a smooth spatial distribution and conserves the same multipole moments. As a result the Fourier coefficients ${~{n}}^{ps}({G}_{s})$ of the pseudo-density for the stars are estimated instead of $~{n}({G}_{s})$ in (41) and then the Coulomb potential part of $~{V}({G}_{s})$ in (6) is written for G_s ≠ 0 by

$$\begin{equation}{~{V}}_{\text{C}}^{ps}({G}_{s})=\frac{8\pi {~{n}}^{ps}({G}_{s})}{{G}_{s}^{2}}.\end{equation} \tag{ 47 }$$

As an origin of energy we choose the average of the interstitial Coulomb potential to ${~{V}}_{\text{C}}^{ps}(0)$.

On the sphere boundaries for the pth atomic group, the Coulomb potentials ${V}_{t}^{p}({a}_{p})$ in (8) are estimated directly from the obtained ${~{V}}_{\text{C}}^{ps}({G}_{s})$ as

$$\begin{equation}{V}_{t}^{p}({a}_{p})=\sum\limits _{{G}_{s}\ne 0}{~{V}}_{\text{C}}^{ps}({G}_{s}){j}_{{L}_{t}}({G}_{s}{a}_{p}){b}_{t}^{p}({G}_{s}).\end{equation} \tag{ 48 }$$

The above equation is derived, if $~{V}({G}_{s})$ changes ${~{V}}_{\text{C}}^{ps}({G}_{s})$ in (6) and then if the lattice-harmonics expansion of the stars,

$$\begin{equation}{\varphi }_{s}(\boldsymbol{r})=\sum\limits _{p=1}^{{n}_{a}}\sum\limits _{t}\;{j}_{{L}_{t}}({G}_{s}r){b}_{t}^{p}({G}_{s}){K}_{t}^{p}(\hat{\boldsymbol{r}}),\end{equation} \tag{ 49 }$$

is used, where L_t is an angular momentum for the tth non-zero lattice harmonics and ${b}_{t}^{p}({G}_{s})$ are the expansion coefficients to be obviously decided from (49).

In the sphere region the Coulomb potentials of V( r ) are solved on the boundary condition of ${V}_{t}^{p}({a}_{p})$ [45] and it can be expressed by

$$\begin{equation}{V}_{\text{C}}(\boldsymbol{r})=\sum\limits _{p}\left(-2{Z}_{p}\left[\frac{1}{r}-\frac{1}{{a}_{p}}\right]\sum\limits _{\nu =1}^{{n}_{p}}\theta ({a}_{p}-{r}_{\nu }^{p})+\sum\limits _{t}{v}_{t}^{p}(r){K}_{t}^{p}(\hat{\boldsymbol{r}})\right),\end{equation} \tag{ 50 }$$

with

$$\begin{align}{v}_{t}^{p}(r)& =\frac{8\pi }{2{L}_{t}+1}\left[\frac{1}{{r}^{{L}_{t}+1}}{\int }_{0}^{r}{n}_{t}^{p}(x){x}^{{L}_{t}+2}\enspace \mathrm{d}x+{r}^{{L}_{t}}{\int }_{r}^{{a}_{p}}{n}_{t}^{p}(x){x}^{1-{L}_{t}}\enspace \mathrm{d}x\right.\\ & \quad -\left.\frac{{r}^{{L}_{t}}}{{a}_{p}^{2{L}_{t}+1}}{q}_{t}^{p}\right]+{V}_{t}^{p}({a}_{p}){\left(\frac{r}{{a}_{p}}\right)}^{{L}_{t}},\end{align} \tag{ 51 }$$

where ${q}_{t}^{p}$ is a multipole moment of L_t for any atom of the pth group,

$$\begin{equation}{q}_{t}^{p}={\int }_{0}^{{a}_{p}}{n}_{t}^{p}(r){r}^{{L}_{t}+2}\enspace \mathrm{d}r.\end{equation} \tag{ 52 }$$

In the last term of the right side in (2), the interstitial XC potentials are calculated from n( r ) of (41) and m( r ) of (42). In order to manage them in practice, beforehand, we make a composite-transformation (CT) method that any function in a real space is transformed by the FFT algorithm to a reciprocal or plane-wave space and subsequently the resulting reciprocal-space functions are symmetrized to ones in a space of the stars. From $~{n}({G}_{s})$ and $~{m}({G}_{s})$, the inverse CT method yields the numerical tables of n( r ) and m( r ) at grid-mesh points of the unit cell and hence the XC potentials are obtained at the same points. Furthermore the XC coefficients ${~{V}}_{\text{XC}}^{ps}({G}_{s})$ for the stars are calculated by the CT method, and the scalar potentials $~{V}({G}_{s})$ in (6) get from $~{V}({G}_{s})={~{V}}_{\text{C}}^{ps}({G}_{s})+{~{V}}_{\text{XC}}({G}_{s})$. The scalar magnetic fields $~{B}({G}_{s})$ in (7), which are oriented along the z axis of the global coordinate, are calculated by the same procedure.

In the sphere region the XC potentials and the magnetic fields are evaluated using n( r ) of (45) and m( r ) of (46) at sampling points inside the spheres and the corresponding ${v}_{t}^{p(\mathrm{X}\mathrm{C})}(r)$ and ${B}_{t}^{p}(r)$ for the lattice harmonics are estimated by the method of least squares so as to be the obtained values in the sampling points [46]. Finally the scalar potentials in the sphere region V( r ) can be determined by adding ${v}_{t}^{p(\mathrm{X}\mathrm{C})}(r)$ into ${v}_{t}^{p}(r)$ in (51).

In this section, we give the details of the DLASTO calculation formulated in the previous section. When one selects a solid to calculate, the input data for the crystal structure is usually set to lattice constants, crystal space group, and Wyckoff parameters including atom types. Using the symmetry operator of the space group, inequivalent atomic sites and their atomic numbers are identified in parallel with checking whether the input data matches the crystal structure, after which the pν atoms are numbered.

To start the DLASTO calculations, the initial crystal densities, which consist of the charge densities n( r ) indicated by (41) and (45) and the magnetizations m( r ) in (42) and (46), are estimated by solving the SPCD equations for (24) and (25). In the sphere region, the SPCD equations are calculated self-consistently under the atomic boundary condition to create the atomic densities for each of the inequivalent atoms. The obtained magnetizations are directly used as the spherically symmetric part ${m}_{0}^{p}(r)$ of m( r ) in (46) and the other parts set zero, because of their atomic-like distributions in solids. While the spherically symmetric part ${n}_{0}^{p}(r)$ of n( r ) in (41) are evaluated from overlapping the obtained charge densities by means of the Löwdin alpha expansion [35] according to the given ${\boldsymbol{\tau }}_{\nu }^{p}$ vectors and the translational vectors in the crystal structure. In the interstitial region, $~{n}(0)$ in (45) sets a constant so as to conserve the total charge in the solid and $~{m}(0)$ in (46) puts a zero constant as an initial value.

After these initial densities have been determined, the potentials in section 2.5 can be calculated. By numerically solving the generalized eigenvalue problem in (35), the energy eigenvalues and eigenvectors expected by these potentials can be obtained, and finally the Bloch function can be obtained. From this Bloch function, new densities can be estimated. The next self-consistent calculation is performed using the new input densities, until the input densities (or the initial densities) and this new output densities agree within a given tolerance. In general, the new input densities are estimated from the previous input and output densities. Here the Broyden's method is often used [47].

The DLASTO method has two cut-off parameters we have to set up in common with the LAPW method. One is an upper limit ℓ_max for the sum of ℓ' in the sphere basis function (21) and another is k_max for the sum of k _i in the interstitial basis function (12). Both of them make a decisive influence on the convergence for basis functions, and eventually for Bloch function. In the APW method, it has been indicated that ℓ_max and k_max are related with values of function j_ℓ (ka_p ) [48], which traces back to the spherical-harmonics expansion of plane waves used to match on the boundaries between the sphere and interstitial regions. In the DLASTO method, it exists explicitly in (32) and (33). As an empirical rule, it is proposed that the d band energies are converged to the order of 0.001 Ry with a criterion k_max a_p = 7 in APW method [48] and in the LAPW method it is a good convergence with k_max a_p ranging from 7 to 9 for all the cases [14]. Although ℓ_max is connected with L_t of the lattice harmonics (A.11) in the full potentials, ℓ_max = 8 is a typical value in the LAPW method [14] or the DLAPW method [22], so as to obey the rule of k_max a_p = ℓ_max in the APW family.

Moreover the DLASTO method has yet another relationship to ζ_p in ${~{\phi }}_{n\ell }^{p}(k,{\zeta }_{p})$ through ${~{\varphi }}_{n\ell s\mu }^{p}({\boldsymbol{k}}_{j},m)$ in (30) and (31). For the 3d orbital, the Fourier component of the radial STO in (17) can be expressed with n = 3, ℓ = 2 and k = k_max by

$$\begin{equation}{~{\phi }}_{32}^{p}({k}_{\mathrm{max}},{\zeta }_{p})=48{k}_{\mathrm{max}}^{2}\frac{{\zeta }_{p}}{{({k}_{\mathrm{max}}^{2}+{\zeta }_{p}^{2})}^{4}}.\end{equation} \tag{ 53 }$$

Thus the computational accuracy of the DLASTO method is partially determined by the choice of the cut-off values of k_max and ζ_p . For the other orbitals, the relation between k_max and ζ_p can be analytically derived from ${~{\phi }}_{n\ell }^{p}({k}_{\mathrm{max}},{\zeta }_{p})=0$ in (17), but in the case of the 3d orbital the optimal solutions of k_max and ζ_p cannot be identified from the function form of (53). Finally the DLASTO method becomes a conditional convergence for k_max and ζ_p .

As an empirical rule, it has been found that in the setting range of k_max a_p = 9–14, the band structure can be well converged irrespective of the value of ζ_p decided from the condition (20) and without a problem of incompleteness of the basis set that frequently occurs in the DLAPW method. It is a larger value than k_max a_p = 8 in the APW family, but it does not increase the computing time, because the Hamiltonian H_NN' and overlapping O_NN' matrix elements in the interstitial region, which are the most time-consuming parts of the DLASTO method, are calculated by the FFT algorithm. Also, by empirically adjusting ζ_p , a smaller value of k_max a_p = 9 is in some cases sufficient to achieve a convergence. As increasing n and ℓ, the convergence tends to occur at small values of k_max a_p due to the property of the function ${~{\phi }}_{n\ell }^{p}(k,{\zeta }_{p})$.

For the spherical basis set (10) we have to choose orbitals {nℓsμ} to every pth atomic group in the solid. In the previous papers [20], it has been suggested that one choice is called 'single ζ' as using a single set composed of orbitals with a different ℓ, and the other is 'double ζ' which is a single ζ plus one more orbitals with an unoccupied states. In the LAPW method, the single and double ζ correspond to a concept of multiple windows. In each of the windows, a distinct set of fixed energies ${E}_{\ell \mu }^{\alpha p}$ is used to handle separate orbitals or semi-core states. Here we have chosen a minimum basis set with single ζ to give an overview of the qualities of the band structure obtained in the DLASTO calculation.

As a test survey, we pick out a typical ferromagnet of body-centered-cubic (bcc) Fe and show the calculated results in the DLASTO method. From here on, we will proceed according to the exchange and correlation potentials of von Barth and Hedin [3]. For the basis set in the sphere region, total 18 orbitals of 3d_2/3, 3d_5/3, 4s_1/2, 4p_1/2 and 4p_3/2 are set up and the other occupied states are calculated as core electrons by solving the SPCD equation or the Dirac equation. As an empirical rule in the 3d transition metals, it has been confirmed that unless the 3p electrons were calculated as valence band electrons, the DLAPW method could not make a spin polarized, as well as in the LAPW method [50]. In the DLASTO method, however, the final result yields the same spin polarization irrespective of whether the 3p states are core electrons or valence electrons. In order to show a wide energy range of the band dispersion, it seems better to add the unoccupied states 4p_1/2 and 4p_3/2, since their bands are distributed near above the Fermi level. The lattice constant 5.33a_B is used in the same values of the previous paper [21] to compare with the DLAPW results, where a_B is the Bohr radius. In the calculation, the value of a_p is set out as 3.3a_B to a certain maximum, and k_max a_p = 12 and ℓ_max = 8 are used. The integration of the densities in section 2.4 is executed by the tetrahedron method at 21 × 21 × 21 grid points in the first octant of the first Brillouin zone, namely at 505 points in the irreducible wedge.

Figure 1 shows the electronic band structure of bcc ferromagnet Fe calculated by the DLASTO method in the full potential scheme (FP-DLASTO), where the Fermi level is indicated in a dashed line. In figure 1, each blue and red color denotes the component intensities of the s and d orbitals, and the rest are represented by black. Thus, the change in color of red and blue reveals the strength of the hybridization between the d and s orbitals. The obtained FP-DLASTO band structure is well agree with the DLAPW one in bcc Fe [21] and also the Fermi surfaces have the same shape as in the DLAPW method. A common feature of the spin-polarized relativistic band structures of the DLASTO and DLAPW methods, other than the spin–orbit splitting, is that since the up- and down-spin states are not in good quantum numbers, a band-splitting occurs when the bands with these two spin states cross each other, which is different from the scalar relativistic band structure, especially in the Fermi-surface shapes. Hence the calculated Fermi surface can explain the angular dependence of the measured dHvA frequencies, though once more details are not given here on the analysis of the Fermi surface [21]. We can see that the FP-DLASTO method works well in the band structure of bcc Fe. If the similar settings were prepared for the 3d transition metals, the FP-LASTO method could determine the band structures within a similar accuracy.

Zoom In Zoom Out Reset image size

Table 1 shows the spin moment μ_S, orbital moment μ_L and SC in bcc Fe in different methods. The spin moment is evaluated from a unit-cell-volume integral of m( r ) in (42) and (46) by

$$\begin{equation}{\mu }_{\text{S}}=\frac{1}{{\Omega}}{\int }_{{\Omega}}m(\boldsymbol{r})\mathrm{d}\boldsymbol{r},\end{equation} \tag{ 54 }$$

and it can be partitioned into the pth atomic group and angular momenta, though the detailed expressions are not given here. It is to be pressed that owing to an orthogonality of plane waves and a spatial symmetry of lattice harmonics, the main contribution of μ_S comes from the spherical term of L_t = 0 in the sphere region and the constant term of G_s = 0 in the interstitial region. Similarly the orbital moment is obtained from the sphere-region basis Ψ_N ( k , r ) in (21) as a projected quantity onto the z-axis of the local coordinate as follows:

$$\begin{equation}{\mu }_{\text{L}}=\sum\limits _{\boldsymbol{k}\lambda }\;{W}^{\lambda }(\boldsymbol{k})\sum\limits _{N{N}^{\prime }}\;{\delta }_{p{p}^{\prime }}{\delta }_{\nu {\nu }^{\prime }}{U}_{N}^{\lambda {\ast}}\left(\frac{1}{{\Omega}}{\int }_{{\Omega}}{{\Psi}}_{N}^{{\dagger}}(\boldsymbol{k},\boldsymbol{r})\beta {\hat{\ell }}_{z}{{\Psi}}_{{N}^{\prime }}(\boldsymbol{k},\boldsymbol{r})\mathrm{d}\boldsymbol{r}\right){U}_{N}^{\lambda },\end{equation} \tag{ 55 }$$

where ${\hat{\ell }}_{z}$ is an operator for orbital angular momentum along the z-axis of the local coordinate. Practically, (55) can be computed using the matrix component of ${\hat{\ell }}_{z}$ as

$$\begin{equation}\langle {\chi }_{s}^{\ell ;\mu }\vert {\hat{\ell }}_{z}\vert {\chi }_{{s}^{\prime }}^{{\ell }^{\prime };{\mu }^{\prime }}\rangle ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{\mu {\mu }^{\prime }}\left(\begin{matrix}\hfill \ell {u}_{\mu }^{\ell }\hfill & \hfill \frac{1}{2}\sqrt{1-{u}_{\mu }^{\ell }}\hfill \\ \hfill \frac{1}{2}\sqrt{1-{u}_{\mu }^{\ell }}\hfill & \hfill (\ell +1){u}_{\mu }^{\ell }\hfill \end{matrix}\right),\end{equation} \tag{ 56 }$$

where the ss' matrix elements for $s={\pm}\frac{1}{2}$ and ${s}^{\prime }={\pm}\frac{1}{2}$ are arranged in the same way as the CG coefficients in (19) using the parameter ${u}_{\mu }^{\ell }$.

For the purpose of investigating a spatial dependence of potentials in the DLASTO method, three different types of potentials: the full (FP-), warped muffin-tin (WM-), and MT potentials are used to calculate the μ_S, μ_L and SC in the same lattice parameters, and the calculated results are shown in the top three lows of table 1. The WM-LASTO calculation uses a spherically-symmetric potential only in the sphere region, and namely it takes into account the same spatial distortion in the interstitial region as in the FP-LASTO one, while the MT-LASTO calculation further makes the interstitial potential constant. Comparing with the three potentials, the spin moments are in agreement within ±0.01 μ_B, while the orbital moment in the FP-DLASTO calculation is slightly larger than in the other ones. It could be seen that the non-spherical potentials in the sphere region are important to enhance the orbital moment in the DLASTO method, though the absolute values are generally small in bcc Fe.

In the four rows or less of table 1, the previous calculated results of the DLAPW method [21] or the other methods [34, 51–54] and the experimental results [55] are shown to compare with the DLASTO results on their calculation precision and trend. In other papers, the relativistic spin polarization theories based on the SPCD equation are often expressed in the abbreviated notation SPR, but note that here we use the same notation, D, for a unified representation. Although there are differences in the exchange potential and correlation potential formulas and lattice constants used, the spin moments have a magnitude of 2.22 μ_B as an average value. The DLASTO results show a spin polarization of about the average value, which is 0.09 μ_B larger than the experimental value. In the orbital moment, it appears that there is almost no change due to the potential shape and the relativistic spin polarization in the KKR and LMTO methods, though they cannot be compared in an exact sense because they were not performed under the same calculation conditions.

The SC of electronic specific heat γ is proportional to the density of states (DOS) at the Fermi level D(E_F) as $\gamma =\frac{1}{3}{\pi }^{2}{k}_{\text{B}}^{2}D({E}_{\text{F}})$, where k_B is the Boltzmann constant. As can be seen from table 1, the SC of the DLASTO method is almost unchanged within the range of 0.08 mJ K⁻² mol⁻¹ with respect to the potential shape, which means that the dispersion of the bands at the Fermi level is unchanged to a corresponding magnitude. Moreover these SC values are about half smaller than the experimental values regardless of the calculation method. The DOS at the Fermi level, or the SC, needs to be improved by a magnitude of about 1 as a dimensionless enhancement factor (=γ_exp/γ_band − 1) in the exchange and correlation potentials or in the fundamental equations, where γ_band and γ_exp are the theoretical and experimental SC, respectively.

In most lanthanide compounds, Fermi surface studies of the dHvA effect have shown that the 4f electrons are localized in the ground state. In short, there are no 4f electrons on the Fermi surface. However the LSDA calculations in any case indicate a heavy-fermion electronic structure with 4f electrons on the Fermi surface: it is crucial to emphasize that the valence state is different from the one indicated by LSDA for the magnetic lanthanide system. Empirically, the angular dependence of the dHvA frequencies observed in magnetic lanthanide compounds can be usually extrapolated, as an alternative, from the topology of the Fermi surface calculated in the La or Y and Sr reference material, not containing 4f electrons, for example, in EuCo₂Si₂ [56]. Among lanthanide magnets, the elemental metal Gd is one of the ferromagnetic materials at room temperature and the only one in the lanthanide series. It crystallizes in the hexagonal-closed-packed (hcp) structure with a Curie temperature of 293.2 K [57]. A general discussion of first principles calculations for Gd metals has been summarized elsewhere [58]. At the onset, an attempt is made to calculate the ferromagnetic state of hcp Gd within LSDA in order to verify the accuracy and trend of the DLASTO method and the obtained results are compared with those of the DLAPW method previously calculated in the same framework [21].

In the LAPW calculation of Gd metal, it has been pointed out that fully occupied 5s and 5p orbitals should be treated as valence electrons including 4f, 5d and 6s orbitals. [50]. According to this criterion, we initially performed the DLASTO calculations with the same configuration of valence electrons as in the previous DLAPW calculations, but we encountered some problems: in the DLASTO method, if the STOs, obtained from ζ_p by (20), are connected to the wave functions of the 5s and 5p orbitals on the sphere, then the smearing contribution to the interstitial region tends to be larger for both orbitals than the one treated as core electrons. The electronic states of Gd metal have actually shown that the lost numbers of these 5s and 5p states from the sphere result in a change in the occupied number of 4f states, eventually giving a large change in the magnetic moment for the Gd ferromagnet. In order to avoid this problem in the DLASTO method, although it is possible to enlarge ζ_p manually to fit the current situation, as an option, it was to set a large radius a_p = 3.2a_B, as allowed by the crystal structure, so that the 5s and 6p orbitals would be spatially well localized inside the sphere as core electrons. In fact, since the 5p electron has only 0.02 electrons in the interstitial region, this is a permissible approximation that can be regarded as core electrons. Therefore we configured a total of 64 STO's in the interstitial region: 4f_5/2, 4f_7/2, 5d_3/2, 5d_5/2, 6s_1/2, 6p_1/2, and 6p_3/2. The basis functions in the sphere region are expanded up to ℓ_max = 8 and the cut-off to the interstitial plane waves is executed by k_max a_p = 14. The lattice constants are set as a = 6.6389a_B and c = 10.6024a_B in the previous results of the DLAPW method [21], which are calculated in terms of the WS radius r_WS = 3.642a_B and c/a = 1.597. The integration of the densities is done by the tetrahedron method at 462 points in the irreducible triangular-prism of the first Brillouin zone of the hexagonal lattice.

Figure 2 shows the electronic structure of the hcp ferromagnet Gd in the DLASTO calculation within LSDA along high-symmetry axis of hexagonal Brillouin zone and the dashed line denotes the Fermi level. The red color in figure 2 illustrates the distribution intensity of the 4f electrons, while the blue color shows the distribution of the 5d electrons and the other contributions are drawn in black. Since the crystal structure of hcp contains two Gd atoms in the hexagonal Bravais lattice, the band structure involves a total of twenty-eight spin-splitting 4f bands, twice the occupation number of fourteen 4f-electrons. As can be seen from the band structure in figure 2, the 4f bands are split above and below the Fermi level into fourteen majority (up-spin) bands and fourteen minority (down-spin) bands with a size of about 4 eV. The magnitude of this splitting is called the magnetic exchange splitting δE_ex, which is of the size of 12 eV in experimental measurements [59, 60]. Although the value of δE_ex in the DLASTO method is very small compared to the experimental value, this band structure is justified by the LSDA calculations: below the Fermi level, there are the nearly dispersionless majority 4f bands at about 4 eV and the large bandwidth s band, and above the Fermi level, there are the minority 4f bands directly above and the widely distributed 5d bands [61]. In particular, for the band structure in the vicinity of the Fermi level, a detailed comparison can be made with the paper reported by Sticht and Kübler [62] using the ASW method including the spin–orbit interaction, and the dispersion of the 5d bands is found to be in very good agreement. The DLASTO method could be concluded to be sufficiently accurate even for the almost dispersionless 4f band if it is calculated with the above indicated cut-off parameters and STO's.

Zoom In Zoom Out Reset image size

Compared with the previous band structure of the DLAPW method [21], the size of δE_ex is the same, but the differences are that majority and minority 4f bands of the DLAPW method have larger dispersion and bandwidth than those of the DLASTO method and in addition the DLAPW Fermi surfaces are formed by the 4f–5d hybridized bands with relatively large dispersion. As seen from figure 2, however, the DLASTO results show that the minority 4f bands are pinned just above the Fermi level, giving an itinerant electron contribution. The occupation number per atom for f orbitals in the sphere, n_f , is shown to be 7.30 by the tetrahedron method. If it is a localized 4f electron with δE_ex = 12 eV as shown in the experiment, then n_f = 7, which means that 0.3 f-electrons are supplied as itinerant electrons in the DLASTO calculations with LSDA. This is the cause of different quantitative values of physical quantities between the DLAPW and DLASTO methods for the Gd ferromagnet, for example magnetic moments, SC and so on.

Therefore, in table 2, the calculated magnetic moments, orbital moments, and SC for the hcp ferromagnet Gd are summarized together with the results from several band theories that can be compared with which the usual first-principles calculations with LSDA have been performed. In the 'method' column of table 2, the initial letter 'D' indicates a spin-polarized relativistic band theory based on the SPCD equation, while the letter '+SO' at the end indicates a relativistic band theory incorporating the spin–orbit interaction into the scalar-relativistic equation by perturbation. In the case of Gd ferromagnet, the magnetic moments are evaluated from the sum of the spin and orbital moments from all orbitals, not just the contribution from the 4f electrons. In the experiment, we only know the definite data that the g-factor is g = 2.00 ± 0.02 [64] and the magnetic moment is 7.63 μ_B [65]. Assuming that the rate of contribution in spin and orbital moments is equal for all orbitals, the magnitude of the orbital moment is estimated to be ±0.08 μ_B. On the other hand, under the condition that the 4f electrons are localized, the orbital moment from the 4f electrons would have disappeared since the seven occupied states are fully spin-polarized, so if we consider only the 0.63 μ_B contribution from non-4f electrons in the same way, the magnetic moment can be estimated to be as close to zero as possible, ±0.006 μ_B. The magnetic moment obtained by the DLASTO method is 7.51 μ_B, which appears to be a reasonable value, but the orbital moment is 0.80 μ_B, which is an extremely large value due to the itinerant contribution from the minority 4f electrons pointed out above. Furthermore, the value of SC is 31.65 mJ K⁻² mol⁻¹ in the heavy-fermion state, which is also a natural conclusion, and does not agree with the experimental results at all.

Although the electronic structure obtained from the with-LSDA-DLASTO calculation cannot be compared with the experimental results, the physical quantities related to the Fermi surface are shown below: the 19th, 20th, 21st, and 22nd bands, counting from the bottom of the valence band, cross the Fermi level and form the Fermi surface. The hole carrier numbers of the 19th and 20th bands are 0.013 and 0.091, and the electron carrier numbers of the 21st and 22nd bands are 0.102 and 0.002. These carrier numbers vary depending on the potential shape and calculation method, but it is clear that the Gd ferromagnet is a compensated metal with the same numbers in the holes and electrons.

Next, let us consider phenomenologically under what conditions we can upgrade the basic equation such as the SPCD equation of the DLASTO method in order to achieve an improvement in the magnitude of δE_ex. The x-ray photoemission spectroscopy of the Gd ferromagnet shows the occupied 4f-state at −7.44 eV from the Fermi level and the unoccupied 4f-state at +4.04 eV [59]. This indicates that for the band structure in figure 2 from the LSAD calculations, the occupied 4f band should be shifted to about 4 eV lower energy and the unoccupied 4f band should be moved to the same value of about 4 eV higher energy. Accordingly it means a localized 4f state in which there are no 4f electrons on the Fermi surface. The basic splitting structure inherent in the SPCD equation needs so as to be converted from a six-to-eight splitting between 4f_5/2 and 4f_7/2 states due to the relativistic spin–orbit interaction to a seven-to-seven splitting between 4f up-spin and 4f down-spin states. In lanthanide trivalent ions, it has been found that the magnetic ground states derived from the SPCD equation in the relativistic jμ representation reproduce the value of Landé g-factor in the LS coupling, where the spin–orbit interaction is negligibly small [26]. This unique feature is due to the Paschen–Back effect, which naturally shifts from the jj-coupling state to the LS-coupling state depending on the magnitude of the internal magnetic field, and is the reason why the DLASTO method can produce reasonable spin-splittings even for Gd ferromagnets, as well as in the DLAPW method.

Considering the spin-splitting structure of the SPCD equation and assuming that an additional constant potential is introduced only for the direct Coulomb term, it is sufficient to introduce a potential of the following form for the SPCD equation of (24) and (25),

$$\begin{equation}{V}^{p}(r)\to {V}^{p}(r)+2\alpha {\delta }_{\ell 3}{(-1)}^{{\delta }_{\mu \frac{7}{2}}}U={V}^{p}(r)+{\Delta}{V}_{\ell \mu }^{\alpha p},\end{equation} \tag{ 57 }$$

depending on the parity 2α = ±1 of 4f_5/2 and 4f_7/2. Here the factor δ_ℓ3 means that this additional potential only acts on the 4f electrons, and the factor ${(-1)}^{{\delta }_{\mu \frac{7}{2}}}$ is introduced according to the Paschen–Back effect. From the above discussion, it is clear that the parameter U in (57) should be set as U = 4 eV ≈ 0.3 Ry. Since the additional potential of ${\Delta}{V}_{\ell \mu }^{\alpha p}$ is a constant value independent of real space, it is equivalent to replacing the band center of the fixed energy ${E}_{\ell \mu }^{\alpha p}$ with ${E}_{\ell \mu }^{\alpha p}-{\Delta}{V}_{\ell \mu }^{\alpha p}$, considering the structure of the SPCD equation between the energy ${E}_{\ell \mu }^{\alpha p}$ and scalar potential V^p (r). Therefore, in the DLASTO method, it is easy to perform the calculations by replacing only the energy term in the Hamiltonian matrix element ${\hat{H}}_{N{N}^{\prime }}^{\text{MT}}$ of (36). Using the explicit expression for ${\hat{H}}_{N{N}^{\prime }}^{\text{MT}}$ described in (A.1) of appendix A, we can substitute ${\varepsilon }_{\ell \mu }^{\alpha p}$ for ${E}_{\ell \mu }^{\alpha p}-{\Delta}{V}_{\ell \mu }^{\alpha p}$. Here we will refer to this relativistic calculation method as LSDA + V. It should be noted here that the LSDA + V calculation is not in a self-consistent procedure, since the additional potential ${\Delta}{V}_{\ell \mu }^{\alpha p}$ in (57) is added only.

Figure 3 shows the electronic structure of hcp ferromagnet Gd calculated by the DLASTO method with LSDA + V along hexagonal symmetry axis of the Brillouin zone. The dashed line denotes the Fermi level and clearly the 4f bands are moving away from the Fermi level, namely being the localized 4f electrons. The Fermi level is crossed mainly by highly dispersive 5d bands. Therefore, the SC value, which is estimated from the DOS of this band structure, is very small, 2.245 mJ K⁻² mol⁻¹, and the enhancement factor λ is very small, 0.64, based on comparison with the experimental value, 3.7 mJ K⁻² mol⁻¹ [66]. Counting from the bottom of the valence band, the 20th and 21st bands form the Fermi surface, which has a hole surface of 13% and an electron surface of 13% of the volume of the Brillouin zone, respectively. This indicates that the localization of the 4f electrons does not change the band index to produce the main Fermi surfaces or the fact that it is the compensated metal.

Zoom In Zoom Out Reset image size

As shown in figure 3, the occupied 4f bands, which have almost flat energies with respect to wavenumber, exist from about 8 eV below the Fermi level in a region of about 2.5 eV width. In the unoccupied state, the flat 4f bands are centered about 4 eV above the Fermi level with about 3 eV width, and nevertheless they form some hybridized and/or mixed bands with the 5d electrons. The occupied 4f bands consist of 14 bands, although the band structure in figure 3 counts only 7 bands, almost flat at this energy scale, because one band in its appearance is actually consisted of two very slightly split states. Likewise, the unoccupied 4f states are made up of 14 bands. From the number of states calculated using the tetrahedron method, it can be confirmed that exactly seven 4f states are occupied per atom. Unfortunately, the magnetization density cannot be determined correctly, because it is not a self-consistent calculation.

Each of the seven apparent 4f bands in figure 3 is identified in terms of the z-component of total angular momentum μ, with the interband intervals widening in increasing energy order for the occupied states and narrowing for the unoccupied states. This splitting structure of the spacing pattern present in the DLASTO calculations is known as Landé interval rule, where the intervals increase in proportion to the magnitude of μ. It is shown that when the internal magnetic field is sufficiently larger than the spin–orbit interaction, this is effectively equivalent to a splitting of the LS coupling with a small spin–orbit interaction. As compared with the experiment [59], however, the calculated 4f interband intervals seem to be quite large. The interband intervals are related to the magnitude of the internal magnetic field in the DLASTO method so that in the practical sense the internal magnetic field allowed for the Gd ferromagnet should be determined in a self-consistent way to balance the large shifts of the 4f bands with LSDA + V. It is important to emphasize that a feasible first-principles reformulation of the underlying equations of relativistic spin polarization, rather than a perturbative approach, requires the introduction of a local potential with both parity and Pachen–Back effect symmetries, as in (57).

Among the actinide compounds, the uranium monopnictides UX (X = N, P, As, Sb and Bi) have a simple NaCl-type crystal structure and their unusual magnetic properties. The space group is Fm3m, and the Bravais lattice is a face-centered cubic (fcc) lattice. Uranium mononitride UN, with the lattice constant a = 4.89 Å, is an antiferromagnet of single- k type I with easy axis of ⟨001⟩ direction and the ordering temperature is 53 K [67]. The SC of the low-temperature specific heat is measured to be 49.6 mJ K⁻² mol⁻¹ [68], and therefore the 5f electrons in UN can be deduced to be itinerant electrons existing on the Fermi surface. First-principles calculations of the electronic structure of UN have been performed using various calculation methods, and detailed discussions of the results have been made elsewhere, including x-ray photoelectron spectroscopy measurements [69]. The application of the DLASTO method to uranium compounds requires caution in several respects. Since the information on band dispersion and Fermi surface of UN has been given by angle-resolved photoemission spectroscopy and the DLAPW calculation results in the paramagnetic state [70], we will restrict our discussion to the paramagnetic band calculation of UN in the actinide system.

Figure 4 shows the charge distributions of the outer orbitals 5f, 6p, 6d and 7s in an isolated uranium atom, where the horizontal axis denotes the distance from the nucleus in units of Bohr radius. These charge distributions are obtained by self-consistently calculating the Dirac equation with LDA, i.e. the equations in (24) and (25) with zero magnetic field, under atomic boundary conditions. The two contributions present in the charge distributions of 5f, 6p and 6d, respectively, illustrate the charge distribution of the two states split by the spin–orbit interaction. That is, for 5f, they are for the 5f_5/2 and 5f_7/2 states. In the case where the radius a_U of the sphere region in the U atom for UN has a value of 3a_B, the vertical solid line in figure 4 corresponds to it. This value of 3a_B is almost the possible maximum value for the UN lattice constant by reducing the radius a_N of the N atom. From this figure, it is clear that the 5f and 6p electrons are well localized and distributed inside the sphere, while the 6d and 7s are spread outside the sphere. Therefore, the 5f and 6p electrons make narrow bands in solids, respectively. On the contrary, the 7s electron is almost unoccupied in solids, and tends to become a band with a very large dispersion, although it is fully occupied in atom. In the DLASTO method, even if ζ_U is fixed at the condition of (20) and a rather large value of k_max is set, it has been confirmed that the band originating from the 7s orbital appears directly in the electronic structure near the Fermi level. Therefore, we extended (20) to include a parameter x as,

$$\begin{equation}{\zeta }_{p}(x)=\frac{n+x}{{a}_{p}},\end{equation} \tag{ 58 }$$

and we investigated the relationship between x and the electronic structure. If x is set to ℓ, we can return to the original value ζ_p (ℓ) = ζ_p . By comparing with the band structure of the DLAPW method calculated previously [70], as a rule of thumb, it is necessary to use an STO with a value of ζ_U around x = 2 or 3 for the 7s electron of U atom.

Zoom In Zoom Out Reset image size

Next, there is a question of whether to approach the 6p_3/2 and 6p_1/2 states of the U atom as core or valence electrons. As can be seen from figure 4, the charge distributions of the 6p_3/2 and 6p_1/2 states are smeared outside the sphere region with a radius of a_U = 3. In fact, the numbers per electron outside the sphere are estimated to be about 0.25 and 0.04 for the 6p_3/2 and 6p_1/2 states, respectively. Therefore, the 6p_3/2 state can only be treated as valence electrons. On the other hand, since the 6p_1/2 state disappears from the sphere with a small amount of electrons, when the calculations are actually done with the 6p_1/2 state treated as core and valence electrons, it is confirmed that there is a small change in the electronic states above and below the Fermi level. In the end, it is concluded that it is better to treat all the 6p states of actinide atoms as valence electrons in the DLASTO method as well as in the DLAPW method. For UN, therefore, the Bloch function is developed using 40 STOs as the basis set in the interstitial region, including 5f_5/2, 5f_7/2, 6p_1/2, 6p_3/2, 6d_3/2, 6d_5/2, and 7s_1/2 in the U atom, and 2s_1/2, 2p_1/2, and 2p_3/2 in the N atom.

Figure 5 shows the band structure of the DLASTO method for UN in the paramagnetic state along the symmetry axis of the fcc Brillouin zone. The settings in the calculation are ℓ_max = 8 and k_max a_p = 10. The radius of the sphere is a_U = 2.8a_B at the U atom, and a_N = 1.7a_B at the N atom. The integration of the densities is done by the tetrahedron method at 280 points in the irreducible wedge of the first Brillouin zone of the fcc lattice. The band structure of the DLASTO method closely reproduces the band dispersion of the DLAPW method, and furthermore it is qualitatively consistent with that of the FPLO minimum-basis calculation, except for a small value of SC [69]. The red and blue colors in figure 5 represent the distribution intensities of U-5f and N-2p electrons, respectively, and the other orbitals are shown in black. The 5f_5/2 and 5f_7/2 states of the U atom exist around the Fermi level, and the U-5f_5/2 bands cross the Fermi level to form the Fermi surface. In addition, some of the U-5f_5/2 states are found to form hybridized bands with the N-2p states. The very flat U-5f_5/2-rich band crossing the Fermi level leads to a large DOS, and the SC of the low-temperature specific heat is estimated to be 18.56 mJ K⁻² mol⁻¹, thus the enhancement factor is relatively small at 1.67. Counting from the bottom of the valence band, the 9th and 10th bands cross the Fermi level and form the hole and electron surfaces of the Fermi surface, respectively. The hole surface occupies 54.5% of the Brillouin zone as holes, and the electron surface occupies 4.5% as electrons. This indicates that UN is an uncompensated metal in the paramagnetic state, and therefore there must be a large Fermi surface with about half the volume of the Brillouin zone.

Zoom In Zoom Out Reset image size

The bands of the occupied states on the symmetry axis Z connecting the symmetry points X and W in the Brillouin zone of the fcc in figure 5 were directly compared with the results measured by angle-resolved photoemission spectroscopy [70]. Although the calculated bandwidth appears to be narrower overall, the number and dispersion of the bands are in very good agreement. In particular, it is noteworthy that one band crosses the Fermi level on the Z axis around the X point. Although the shape of the Fermi surface is not shown here, it is very similar to that of the DLAPW method, and thus the shape of the Fermi surface around the X point is consistent with the Fermi surface mapping of angle-resolved photoemission.