Relativistic full-potential multiple scattering theory: An ab initio method and its applications

TL;DR

This paper introduces a full-potential relativistic KKR method that directly solves Dirac equations for extit{ab initio} calculations, eliminating charge density pathologies and enabling accurate studies of complex materials.

Contribution

It presents a novel full-potential implementation of the relativistic KKR method that improves accuracy for non-spherical potentials in extit{ab initio} calculations.

Findings

Successfully applied to polonium crystal structures and properties

Accurately modeled noble metals and relativistic effects

Eliminated charge density pathologies at small radii

Abstract

The Green function plays an essential role in the Kohn-Korringa-Rostocker (KKR) multiple scattering method. In practice, it is constructed from the regular and irregular solutions of the local Kohn-Sham equation and robust methods exist for spherical potentials. However, when applied to a non-spherical potential, numerical errors from the irregular solutions give rise to pathological behaviors of the charge density at small radius. Here we present a full-potential implementation of the relativistic KKR method to perform \textit{ab initio} self-consistent calculation by directly solving the Dirac differential equations. The pathology around the origin is completely eliminated with the help of an efficient pole-searching technique. This method is utilized to investigate the crystal structures of polonium and their bulk properties. The noble metals are also calculated, both as a test of our method and to study the relativistic effects.

Figures

Tables

Table 1. Table 1: The poles of the single-site Green function of Po

	real part	imaginary part
$s_{1/2}$	-0.68202472D+00	0.16325374D-15
	-0.68202472D+00	-0.21372993D-15
$p_{1/2}$	-0.35205135D-01	-0.12391130D-13
	-0.35205135D-01	0.11887440D-13
$p_{3/2}$	0.66052385D-01	-0.37929455D-01
	0.66052385D-01	-0.37929455D-01
	0.66052385D-01	-0.37929455D-01
	0.66052385D-01	-0.37929455D-01

Table 2. Table 2: Comparison of the lattice constant ( a 𝑎 a ) and bulk modulus ( B 0 subscript 𝐵 0 B_{0} ) calculated with the experimental and theoretical data in the literature.

	a(Å)	$B_0$ (GPa)	Ref.
Exp.	3.345(2)		[33]
	3.359(1)		[34]
LDA, PP	3.28	56	[35]
GGA+SO, FLAPW	3.34		[32]
LDA+SO, FLAPW	3.334	42.3	[37]
LDA+SO, FLAPW	3.323	47.35	[41]
TB+SO	3.29	51	[40]
LDA+Dirac, MST	3.327	55.1	This work

Table 3. Table 3: Comparison of calculated lattice constants and bulk moduli of noble metals with the experimental values. The LDA is employed in all calculations. Fully relativistic schemes are used except for the last two columns. In column Au NR subscript Au NR \rm{Au_{NR}} all electrons are calculated with nonrelativistic schemes, and relativistic schemes are utilized only for valence electrons in column Au R − Core subscript Au R Core \rm{Au_{R-{Core}}} .

	Cu	Ag	Au	${\mathrm{Au}}_{\mathrm{NR}}$	${\mathrm{Au}}_{R-\mathrm{Core}}$
Lattice constants (a.u.)
This work	6.65	7.55	7.60	8.03	8.17
Experiment	6.84	7.72	7.71
Bulk Modulus (GPa)
This work	191	141	208	112	124
Experiment	137	101	180

Equations

$$\rho (\mathbf{r})=-\frac{1}{\pi }\mathrm{Im}\mathrm{Tr}{\int }_{E_b}^{E_F}G(E,\mathbf{r},\mathbf{r})dE,$$

$$N(E)={\int }_{E_b}^{E}n(E^{\prime })d{E}^{\prime },$$

$$n(E)=-\frac{1}{\pi }\mathrm{Im}\mathrm{Tr}{\int }_{\Omega }G(E,\mathbf{r},\mathbf{r})d\mathbf{r}.$$

$$G(E,\mathbf{r},{\mathbf{r}}^{\prime })=$$

$$G_s(E,\mathbf{r},{\mathbf{r}}^{\prime })=$$

$$G_m(E,\mathbf{r},{\mathbf{r}}^{\prime })=$$

$$G_s(E,\mathbf{r},\mathbf{r}{)}_{|E\approx E_n}\approx \frac{{\psi }_n(E,\mathbf{r}){\psi }_n^{†}(E,\mathbf{r})}{E-(E_n-i{\lambda }_n)},$$

$${\int }_{\Omega }{\psi }_n(\mathbf{r}){\psi }_n^{†}(\mathbf{r})d\mathbf{r}=1,$$

$$n(E{)}_{|E\approx E_n}\approx$$

$$F(E)=$$

$$\underline{S}(E)=\left[-i\underline{s}(E)-\underline{c}(E)\right]{\left[i\underline{s}(E)-\underline{c}(E)\right]}^{-1},$$

$$\underline{J}(E)=i\underline{s}(E)-\underline{c}(E).$$

$$\det \left[\underline{J}(z_p)\right]=0.$$

$$\underline{J}(z)=\underline{J}({ϵ}_0+\lambda )={\underline{J}}^{(0)}({ϵ}_0)+{\underline{J}}^{(1)}({ϵ}_0)\lambda +{\underline{J}}^{(2)}({ϵ}_0){\lambda }^2,$$

$$\underline{A}(\lambda )={\left[J^{(2)}({ϵ}_0)\right]}^{-1}\underline{J}({ϵ}_0+\lambda ).$$

$$\underline{A}(\lambda )={\lambda }^2-\underline{B}\lambda +\underline{C},$$

$$\underline{B}=-{\left[{\underline{J}}^{(2)}({ϵ}_0)\right]}^{-1}{\underline{J}}^{(1)}(ϵ);\underline{C}={\left[{\underline{J}}^{(2)}({ϵ}_0)\right]}^{-1}{\underline{J}}^{(0)}(ϵ),$$

$$\underline{A}(\lambda )={\lambda }^2-\underline{B}(\lambda -\underline{D}),$$

$$\underline{D}={\underline{B}}^{-1}\underline{C}=-{\left[{\underline{J}}^{(1)}({ϵ}_0)\right]}^{-1}{\underline{J}}^{(0)}(ϵ).$$

$$\det \left[\underline{A}(\lambda )\right]=\det \left[{\lambda }^2-\underline{B}(\lambda -\underline{D})\right]=0,$$

$$\det \left[\begin{array}{cc}\lambda -\underline{D}& -\underline{I}\\ {\underline{D}}^2& \lambda -\underline{B}+\underline{D}\end{array}\right]=0.$$

$$\left[\begin{array}{cc}\underline{D}& \underline{I}\\ -{\underline{D}}^2& \underline{B}-\underline{D}\end{array}\right].$$

Full text

1 Introduction

Multiple-scattering theory (MST) underpins a number of widely used methods for solving the electronic structure problem in periodic solids, all of which have their origins in the KKR method originally introduced by Korringa[1] in 1947 and independently re-derived by Kohn and Rostoker[2] in 1953. Two features of MST distinguish it from conventional Raleigh-Ritz variational approaches. Firstly, it naturally yields a separation between the single-site potential scattering and structural arrangement (positions) of the individual scatterers. Secondly, in the framework of Density functional theory (DFT), it provides an explicit expression for the Green function of the system, which can then be used to calculate the charge and spin densities without explicit calculation of the wavefunctions that are the focus of Raleigh-Ritz methods. The availability of the Green function makes MST a versatile tool that can be easily combined with other methods to investigate more complex systems. For example, by applying the Dyson’s series expansion to the Green function, defects and impurities in an otherwise perfect crystal can be investigated [3]. Another example is the KKR-CPA method, which is based on a combination of MST with the coherent potential approximation (CPA) [4, 5, 6] to calculate the configurationally averaged properties of disordered systems, such as random alloys [7, 8] and the disordered local moment state of metallic magnets [9]. A more recent development is in studying of strongly correlated systems, where the MST Green function can be readily used in conjunction with the GW approximation[10] or the dynamical mean field theory (DMFT) [11]. Moreover, the real space formulation of MST [12] has demonstrated essentially ideal linear scalability on current supercomputing architecture[13], and, as a result, can be employed to study solid state systems with tens of thousands of atoms.

The originally formulated MST solved the Schrödinger equation within the muffin-tin (MT) potential approximation, where the potential is assumed to be spherically symmetric within the muffin-tin spheres and constant in the interstitial region [1, 2]. While the muffin-tin approximation generally works well for systems dominated by metallic bonding, it cannot properly describe a wide range of systems where the asymmetries of the effective potential [14] play an important role, such as surfaces, two-dimensional materials, and systems with directional covalent bonding. In addition, because the Schrödinger equation is nonrelativistic, it cannot properly describe systems where relativistic effects are important. In particular, it doesn’t account for the spin-orbit coupling (SOC), a subject currently of great interest due to its role in many technically important phenomen, such as magnetocrystalline anisotropy, Rashba effect, and magnetic Skyrmions[15]. To take into account the relativistic effects, a common practice is to treat the relativistic kinematic effects with the scalar-relativistic approximation [16], and include the SOC in a perturbative second-variational way. However, this strategy could be problematic for heavy-element systems where SOC is not small compared crystal field splitting. To take into account both relativity and the full shape dependence of the crystal potential on an equal footing, the original MST formulation must be extended to a full-potential, Dirac equation, based theory, and indeed much work has been done in this regard by a number of groups[17, 18, 19, 20, 21, 22, 23].

In MST the Green function is constructed from the regular and irregular solutions of the Kohn-Sham equations. In contrast to the MT scheme, a persistent problem in standard implementation of full-potential MST is that numerical errors in the irregular solutions are very difficult to control near the origin[24]. As a result, the charge density calculated from the Green function exhibits pathologies which can extend to a sizable fraction of the muffin-tin radius. The practice employed by Huhne et al [20] is to drop the non-spherical components of the potential within a cutoff radius $r_{ns}$, with the argument that close to the nucleus the non-spherical contributions to the potential are small. However, for an accurate determination of the Hellmann-Feynman forces, these non-spherical components of the potential are important and cannot simply be discarded. In Ref. [25], a modified single-site Green function is proposed to avoid directly using the irregular solutions, which, however, still show up a volume integration in the expression of the Green function. In Ref. [26] it is proposed to use a sub-interval technique to systematically decrease the numerical error by reducing the step size when approaching the origin. This method still requires spherical potential approximation within a small radius, and is also less effective for larger $l_{max}$, which signifies the angular momentum cutoff of the solutions.

To completely get rid of this pathology, we use an approach developed from the work of Rusanu et. al.[24], where the use of irregular solution is avoided by splitting the full Green function into two parts and carrying out the energy integration along different contours on the complex plane. One major difference between our method and the one in [24] is in the energy integration of the single-site Green function, where the presence of the single-site resonance states and shallow bound states make a naive integration scheme laborious on the real energy axis. By making use of an efficient pole-searching algorithm of the authors, Y. Wang, we find this real axis energy integral can be accomplished much faster. Furthermore, because it explicitly identifies both bound and virtual bound electron states, our method provides an excellent framework for implementing schemes, such as LDA+U [27] and self-interaction correction (SIC) [28, 29], aimed at correcting local approximations to the DFT for the effects of strong correlation.

In the following section we will first explain the scheme that is used for performing contour integration of the MST Green function, then we will show how the poles of the single-site Green function can be used to facilitate the energy integration of the shallow bound states and the resonance states. Details of this pole-searching technique are presented in the appendix. In section 3, polonium is used as an example to demonstrate our method. The lattice constants, bulk modulus and crystal structures of Po are studied and compared with results from other methods. In section 4, the density of states and bulk properties of copper, silver, and gold are calculated as a further test of our method and to quantify the increasing impact of relativistic effects.

2 Methods

The two physical quantities of most interest in the present context are the integrated density of states $N(E)$ and the charge density $\rho(\mathbf{r})$. In a typical ab initio DFT calculation, these quantities need to be evaluated at each self-consistent loop to determine the new Fermi energy and effective potential. In MST, the charge density is obtained from the energy integral of the Green function,

[TABLE]

where $E_{b}$ is the bottom of the conduction band, $E_{F}$ is the Fermi Energy. The integrated density of states (IDOS) is given by the energy integral of the density of states (DOS) $n(E)$

[TABLE]

and $n(E)$ is calculated from the volume integral of the Green function

[TABLE]

The Green function is obtained by solving the Dirac-Kohn-Sham equation and the details have been given in Ref. [23]. Note that the preceding equations involve an energy integration along the real energy axis in order to obtain $N(E)$ and $\rho(\mathbf{r})$. Unfortunately, for bulk materials a simple energy integration on the real axis turns out to be infeasible due to the dense set of poles in the corresponding multiple scattering Green function. One resolution of this problem is to carry out the integral along a contour in the complex energy plane [30], with the observation that the Green function is holomorphic except for poles at the bound states and a cut on the real axis starting at $E_{b}$. Because the DOS becomes increasingly smooth the further the contour is distorted into the complex plane, this method has been found to be very efficient. Indeed, deploying Gaussian quadrature integration method, only a few dozen energy points are needed to reach a high accuracy. In practice, however, implementation is hindered by the presence of the irregular solutions in the expression of the Green function, which is commonly written as [23]

[TABLE]

where $\Lambda$ stands for the pair of relativistic angular-momentum indices ($\kappa$,$\mu$), $Z_{\Lambda}(E,\mathbf{r})$ and $Z^{+}_{\Lambda}(E,\mathbf{r})$ are the right-hand and left-hand regular solutions, respectively, $\mathcal{J}^{+}_{\Lambda}(E,\mathbf{r}^{\prime})$ are the left-hand irregular solutions, and $\tau_{\Lambda\Lambda^{\prime}}(E)$ are the scattering-path matrices. The irregular solutions are singular at the origin and are obtained by integrating inward from outside of the bounding sphere. Unfortunately, using standard numerical integration algorithms, the irregular solutions typically have unacceptable numerical errors near the origin, which then results in the aforementioned pathology in the charge density. . The reason for this instability is as follows: In full-potential scheme, the non-spherical potential components couple solutions of the differential equation having different $l$ indices. Near the origin, the irregular solutions diverge as $r^{-l-1}$, and the coupling of this divergence to that of channels of higher $l$ amplifies the numerical round-off error in the irregular solutions; an effect that is further amplified as the $l_{max}$ use in the differential equation solver is increased.

In our approach, the elimination of this pathological behavior of the Green function is accomplished by spitting the Green function into the single-site scattering part $G_{s}$ and the remaining multiple scattering part $G_{m}$, as suggested in Ref. [24],

[TABLE]

The explicit expression of the $G_{s}(E,\mathbf{r},\mathbf{r}^{\prime})$ is

[TABLE]

where $t_{\Lambda\Lambda^{\prime}}(E)$ is the single-site $t$ matrix. For simplicity of the discussion, we only consider the case of one atom per unit cell. The explicit expression of $G_{m}(E,\mathbf{r},\mathbf{r}^{\prime})$ is

[TABLE]

Note that the $Z\mathcal{J}$ term in $G_{s}(E,\mathbf{r},\mathbf{r}^{\prime})$ is real for real energies. Because the single-site DOS and charge density involve only the imaginary part of $G_{s}(E,\mathbf{r},\mathbf{r}^{\prime})$ this term can be ignored. As a result, carrying out the required integration over real energy axis obviates the need to evaluate irregular solutions. As for the multiple scattering contribution, because both $G_{s}(E,\mathbf{r},\mathbf{r}^{\prime})$ and $G(E,\mathbf{r},\mathbf{r}^{\prime})$ are holomorphic on the upper half-plane, it can be easily integrated along a semi-circle contour, as shown in Fig. 1.

To efficiently evaluate the energy integral of $G_{s}(E,\mathbf{r},\mathbf{r}^{\prime})$, two obstacles must be overcome. The first one is to properly account for the contribution from sharp resonance states, examples of which are the $d$-state resonances in noble metals [23] and the $p$-state resonances of polonium, as shown in Fig. 3. In relativistic full potential schemes, these resonance peaks usually get sharper and split into multiple peaks due to spin-orbit coupling and crystal field splitting, and make a straightforward energy integration even more prohibitive. The second difficulty is to carry out the energy integration of the shallow bound states, which have poles right on real axis and make a direct numerical integration unfeasible. These shallow bound states show up, for example, in our calculation of polonium, as listed in Table 1.

The resonance peaks originate from the poles located at the forth quadrant in the complex plane. As will be explained in the following, these peaks are well approximated by the Lorentzians, and the energy integration on the positive axis can be carried out efficiently with a weighted sampling technique. To accomplish this, the single-site Green function in the neighborhood of a scattering resonance is first approximated as

[TABLE]

where the complex resonance energy has be written in terms of the real and imaginary part as $E_{n}-i\;\lambda_{n}$. By substitution of equation (8) into equation (3) and using the normalization condition of the wave functions

[TABLE]

then the density of states around $E_{n}$ becomes

[TABLE]

which is exactly a Lorentzian function. The values of $E_{n}$ and $\lambda_{n}$ are determined using the pole-searching technique detailed in section A. Now that the approximate behavior of the resonance peaks are known, we can construct a weighted energy mesh to carry out the integration, i.e., an energy mesh that is densest at the resonance peaks. To use this method, we need to find an appropriate cumulative distribution function $F(E)$. Here it is chosen to be

[TABLE]

where the first part is the integral of the Lorentzian function and the second part is to account for the non-resonance (free-electron) states, with V being the volume of the unit cell. The weighted energy mesh is obtained by uniformly choosing points $F_{i}$ between $F(0)$ and $F(E_{F})$, then solving the inverse of $F(E)$ with, for example, bisection method.

For shallow bound states, the above weight sampling technique no longer works because what was a Lorentzian function at positive energy evolves into a Dirac delta function at negative energies. We instead carry out the energy integration on small contours that encircle the negative energy poles ( see Fig. 1). The reason that radius of the contour must be small is to reduce the error caused by not carrying out the integration strictly on real axis. In our experience, this method yields accurate results when the radius of contour chosen to be $10^{-4}$ Ry.

The numerical parameters used in the calculations are as follows: The angular momentum cut-off of the wave function is chosen to be $l_{max}=4$ and 4096 k-points are used in the first Brillouin zone unless otherwise specified. The number of energy points used in the weighted sampling integration is 100. Thirty Gaussian energy points distributed on a large semi-circular contour that encompasses the full valence band are used to integrate $G_{m}(E,\mathbf{r},\mathbf{r}^{\prime})$ and 5–10 energy Gaussian points are used for the small contours around shallow bound states.

3 Polonium

Polonium is an element that is extremely toxic and highly radioactive, and experimental data on its physical properties is scarce. One distinctive property of Po is that it is the only element to crystallize in the simple cubic (sc) structure at ambient conditions. The simple cubic structure has a low atomic packing factor, and is generally considered unstable, both from the point of views of elastic stability [31] and Peierls instability [32]. In addition to the simple cubic phase Po ($\alpha$-Po), a second, slightly distorted, rhombohedral phase, Po ($\beta-$Po), is also found to exist at elevated temperatures. Both the two crystal structures of Po are studied by Beamer and Maxwell[33] using X-ray diffraction, and the lattice constants are reported to be $a=3.345(2)\,\mathrm{\AA}$ for $\alpha$-Po and $a=3.359(1)\,\mathrm{\AA}$ for $\beta$-Po. Later Sando and Lange[34] redetermined the lattice constants using purer Po sample, with $a=3.359(1)\,\mathrm{\AA}$ for $\alpha$-Po and $a=3.368(1)\,\mathrm{\AA}$ for $\beta$-Po.

With the development of electronic structure calculation methods, a number of theoretical studies have been carried out to explain why Po has a stable sc structure. There are still debates on this question, but the general consensus is that relativistic effects play an important role. Using the pseudopotential (PP) method, Kraig et al. [35] showed that sc structure has the lowest total energy among a number of lattice configurations, including face-centered cubic structure (fcc) and body-centered cubic structure (bcc). Min et al. [32, 36] ultilized the full-potential linearized augmented plane wave (FLAPW) method implemented in the WIEN2K package and found the sc structure is due to SOC. Also using FLAPW, Legut et al. [37, 38] investigated the total energy profile of Po and analyzed contributions from different relativistic terms, and concluded instead that it is the relativistic mass-velocity and Darwin terms that stabilize the sc structure. This controversy is mainly due to the tiny total energy difference (especially for all-electron calculation) between the sc structure and trigonal structure of Po, which is at the order of 0.1 mRy, and can be affected by factors such as accuracy of exchange-correlation functional, convergence of the total energy with respect to numerical parameters (for example, we found the convergence with respect to radial grid is important because for different structures the grid generated can be different.), or even which experimental value of the lattice constant to use in the band and phonon dispersion calculations [39]. Compared to the second-variational implementation of SOC in WIEN2K, the fact that our method includes the SOC intrinsically provides some advantage, but in this paper we do not intend to address the question on which relativistic effects contribute more to the stabilization of sc-Po, mainly because of the aforementioned tiny total energy difference. Instead, we use Po as an example to demonstrate our method, and compare the calculated physical quantities with experiment and other calculation results.

The electron configuration of Po is [Xe]$4f^{14}5d^{10}6s^{2}6p^{4}$. The core and semi-core electrons are calculated by solving the Dirac equation for the spherical (l=0) component of the potential. The valence electrons, i.e., $6s$ and $6p$ states, are calculated with the FP-MST method of this paper. As described in section 2, the Green function is split into a single-site part and a multiple scattering part. The poles of the single-site Green function are computed and listed in Table 1, and the indices in the first column signify the corresponding predominate angular momentum quantum numbers of each pole. The first two $s$ poles and the following two $p$ poles have negative real part and negligible imaginary part, which are the characteristics of bound states. The other four $p$ poles are in the forth quadrant of the complex plane and correspond to resonance states. The splitting of the six $p$ electrons into a doubly degenerate and weakly bound $p_{1/2}$ core state and a $p_{3/2}$ fourfold degenerate resonance state just above the energy zero is due to SOC. Note that the cubic symmetry in the potential does not break the degeneracy of either of these spin-orbit split manifolds. The spherical components of the charge densities of the $s_{1/2}$ and $p_{1/2}$ bound states are shown in Fig. 2. Clearly, the charge densities demonstrate correct number of nodes. Note that these nodes result from the orthogonality of the valence states to the core states and occur for a radius as small as 0.1 a.u.. As in Ref. [24], this is a further example where the previous method of simple polynomial extrapolation the charge density, over some significant fraction of the muffin-tin radius, will miss the effects of these undulations.

The single-site DOS and the total DOS of the system are shown in Fig. 3; the latter is in good agreement with the scalar relativistic + SO interaction calculation result of Ref. [37]. The total energies and the fitted equation of state Po with sc, bcc and fcc structures are shown in Fig. 4. The lattice constant for the sc phase is found to be $a=3.327$ Å, which, as can be seen in Table 2, is in good agreement with the experimental value $a=3.359$ Å. The bulk modulus is determined to be $B_{0}=55.1$ GPa, but unfortunately no reliable experimental value of $B_{0}$ is known. In Ref. [40], it is claimed that the experimental bulk modulus of $\alpha$-Po is 26 GPa, but no reference is given. In Ref. [41] a more complete list of the theoretical results of $\alpha$-Po is given, from which we see all bulk moduli calculated theoretically are in the range of $40-60$ GPa, which is much larger than the 26 GPa quotation. We therefore conclude that the bulk modulus found here is reasonable and consistent with the values in the literature.

4 Noble metals

In contrast to polonium, the noble metals have been thoroughly investigated both theoretically and experimentally. In this section, we use them as examples to further test our method. The lattice constants and bulk moduli calculated here are shown in table 3 together with experimentally measured values. The results are in excellent agreement with the ones calculated with conventional KKR method [42]. Compared to the experimental data, we underestimate the lattice constants and overestimate the bulk modulus, which is a well-known characteristic of LDA. As demonstrated in Ref. [42], by employing GGA instead of the LDA functional, much better agreement with the experimental data can be obtained for noble metals.

To demonstrate the impact of relativity, we calculated the bulk properties of Au with two other schemes. In the nonrelativistic (NR) scheme all relativistic effects are ignored, and in the relativistic core (R-Core) scheme relativistic effects are included only for the core and semi-core electrons. The lattice constants and bulk modulus are shown in the last two columns of Table 3 and the total energy vs lattice constants plot is shown in Fig. 5. From these results, it is clear that relativistic effects need to be included for all the electrons, core and valence, when calculating the ground state properties of elements as heavy as Au.

The impact of relativity is also evident in the DOS plots of Fig. 6. As a point of reference, the fully relativistic DOS are in good agreement with the results in Ref. [43] for the same set of systems. As can be seen, in all nonrelativistic calculations, the DOS have five main peaks, which is a result of the crystal field symmetry. Under the increasing influence of spin-orbit coupling, these peaks further split into multiple sub-peaks. Overall, the differences between relativistic and nonrelativistic results also grow as the atomic number increase. For silver and gold, the $d$-bands obtained in the relativistic calculation broaden significantly with the result that the top of the $d$-band is much closer to the Fermi energy than in Cu. This effect is the result of the relativistic contraction of the inner shell $s$ electrons, whose relativistic effects are more significant than the $d$-electrons. Actually, this is the well-known explanation of the color of gold; namely, that relativistic effects, decrease the transition energy between $5d$ and $6s$ states which then absorb blue light making the reflected light appear golden to us. The relativistic effects in silver are small compared to gold, therefore Ag still reflects all the visible wavelengths and appears “silver”. Concerning the energy difference between top of the $d-$band and the Fermi energy in the relativistic plots, we find they are small compared to those obtained from photoemission measurements[44, 45]. Again, this is a typical feature of DFT calculations, and is shown to be largely corrected [46] by using GW method to account for self-energy effects.

5 Conclusions

We have implemented a new approach to full-potential relativistic MST. By splitting the Green function into two parts and carrying out the energy integration along different contours, this method requires no evaluation of the irregular solutions and therefore is free of any pathology of the charge density near the origin. By explicitly searching the poles of the single-site Green function, we devised an efficient integration scheme to solve the numerical problems caused by the shallow bound states and the resonance states. The density of states and bulk modulus of polonium are calculated, with the lattice constant found to be $a=3.327$ Å, and the bulk modulus $B=55.1$ GPa, which yield excellent agreement with experimental data and results from other methods. As a test of our code, we also calculated the DOS and bulk modulus of Cu, Ag, and Au, and discussed the impact of relativistic effects.

6 Acknowledgements

This work was sponsored by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Material Sciences and Engineering Division. This research used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575. Specifically, it used the Bridges system, which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC). This research also used the resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy.

Appendix A Pole searching technique

The pole searching technique is similar to the one used in the quadratic KKR (QKKR) method to quickly finding the poles of the KKR matrix [47] corresponding to the energy eigenvalues of the electronic band structure. In scattering theory, the bound states and resonance states correspond to the poles of the S-matrix $\underline{S}(E)$ on complex energy plane, which can be written as

[TABLE]

where $\underline{s}(E)$ and $\underline{c}(E)$ are the sine and cosine matrices[23]. To find the poles of S-matrix we only need to identify the zeros of the Jost matrix $\underline{J}(E)$, which is given by

[TABLE]

In scattering theory, the Jost matrix is actually a more fundamental quantity than the S matrix because it has no redundant zeros [48]. To efficiently determine the zeros of the Jost matrix, a linear algebra method is used. Let us consider a square matrix $\underline{J}(z)$ of size $L\times L$ (in relativistic case, $L=2(l_{max}+1)^{2}$), and we need to find its zeros, $z_{p}$, such that

[TABLE]

For the present purposes we are only interested in the poles, $\epsilon_{p}$, that are close to the real energy axis corresponding to either resonance states $(\epsilon_{p}>0)$ in the valence band or the shallow bound states $(\epsilon_{p}<0)$. Accordingly, we choose an energy, in the neighborhood of a pole and perform a quadratic expansion of the matrix around $\epsilon_{0}$ as follows,

[TABLE]

where $\lambda=z-\epsilon_{0}$. In analogy with the terminology used in the quadratic KKR method, we refer to $\epsilon_{0}$ as the pivot energy for the expansion of $\underline{J}(z)$. To find the zeros of this quadratic equation, we consider an alternative matrix,

[TABLE]

By multiplying $\left[\underline{J}^{(2)}(\epsilon_{0})\right]^{-1}$ on both sides of equation (15), we get the following expansion,

[TABLE]

where

[TABLE]

Obviously, $\underline{A}(\lambda)$ has zeros at the same energies, $\lambda_{p}=\epsilon_{p}-\epsilon_{0}$, as does $\underline{J}(\epsilon_{0}+\lambda)$. We now rewrite equation (17) as follows

[TABLE]

where

[TABLE]

The zeros of $\underline{A}(\lambda)$ can now be found using the determinantal equation

[TABLE]

or, equivalently,

[TABLE]

The eigenvalues $\lambda_{p}$ (p = 1, 2,…2L) that satisfy this secular equation can be quickly found by diagonalizing the following matrix:

[TABLE]

The zeros of matrix S(z) are thus at $z_{p}=\lambda_{p}+\epsilon_{0}$. The zeros of matrix S(z) are thus at $z_{p}=\lambda_{p}+\epsilon_{0}$. In practice, a simple way of computing the quadratic expansion coefficient matrices in equation (15) is to calculate $\underline{J}(z)$ at three energy values: $\epsilon_{0}-\lambda,\epsilon_{0},\epsilon+\lambda$, with $\lambda$ an small value, then solve the quadratic expansion equations. To search the poles over the full interval $E_{b}$ and $E_{F}$, a set of panels on the real energy axis is set up and $\epsilon_{0}$ is chosen at the center of each panel to obtain a first approximation to $\epsilon_{p}$. The accuracy of the pole location can be systematically improved through iteration by progressively decreasing the energy window around the pivot energy used to set up the pole location eigenvalue equation.