Benchmark of density functional theory for superconductors in elemental materials
BBenchmark of density functional theory for superconductors in elemental materials
Mitsuaki Kawamura, ∗ Yuma Hizume, and Taisuke Ozaki Institute for Solid State Physics, The University of Tokyo, Kashiwa 277-8581, Japan Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan (Dated: March 4, 2020)Systematic benchmark calculations for elemental bulks are presented to validate the accuracyof density functional theory for superconductors. We developed a method to treat the spin-orbitinteraction (SOI) together with the spin fluctuation (SF) and examine their effect on the supercon-ducting transition temperature. We found the following results from the benchmark calculations:(1) The calculations, including SOI and SF, reproduce the experimental superconducting transitiontemperature ( T c ) quantitatively. (2) The effect by SOI is small excepting a few elements such asPb, Tl, and Re. (3) SF reduces T c s, especially for the transition metals, while this reduction is tooweak to reproduce the T c s of Zn and Cd. (4) We reproduced the absence of superconductivity foralkaline (earth) and noble metals. These calculations confirm that our method can be applied to awide range of materials and implies a direction for the further improvement of the methodology. I. INTRODUCTION
The first-principles calculation of the superconduct-ing properties such as the transition temperature ( T c )and the gap function is of great interest to explore newmaterials as well as to understand the physical mecha-nism of known superconductors. Density functional the-ory for superconductors (SCDFT) is one of the frame-works for such calculations; this method enables us toperform fully non-empirical simulations in the supercon-ducting phases at a reasonable computational cost. Theanisotropic Migdal-Eliashberg (ME) equations and theMcMillan’s formula which is the parametrization of thesolution of the ME equations can also be used to estimate T c . However, to solve the Migdal-Eliashberg equations ,we need to perform the summation of the Matsubara fre-quencies and this summation requires a substantial com-putational cost. Since the McMillan’s formula involvesan adjustable parameter to evaluate the effect of theCoulomb repulsion, the formula cannot compare the T c sof a wide range of materials. In SCDFT, we can treat theelectron-phonon interaction, the electron-electron repul-sion, and the spin-fluctuation (SF)-mediated interaction in a first-principles manner. SCDFT has been appliedto various kinds of materials such as elemental materials(Al, Nb, Mo, Ta, Pb) , MgB , graphite intercalations ,Li under high pressure , H molecule solid , hydro-gen compounds , and FeSe . On the other hand, themethodological improvements have also been proposedto include the anisotropic electron-phonon interaction,plasmons , spin-fluctuation , and the spin-orbit inter-action (SOI) .However, the accuracy of the current approximatedfunctional of SCDFT and the effects of SOI and SFhave not been verified systematically although such ver-ification is highly desirable before applying this methodto a wide range of materials. Such a high-throughputcalculation was performed, for example, in the ex-ploration of low-thermal-conductivity compounds usingfirst-principles calculations together with the materials informatics . A benchmark is also a useful tool usedto find a guideline for improving the theory and ap-proximations of the superconducting density functional.For this purpose, in this paper, we are presenting thebenchmark calculations of SCDFT. As benchmark tar-gets, we have chosen the simplest superconducting andnon-superconducting materials, i.e., elemental materials;each material in this group comprises a single element.The particular computational cost is relatively low be-cause most materials in this group contain only one ortwo atoms in the unit cell. Moreover, we can see the ef-fects of the chemical difference and the strength of theSOI of each element.This paper is organized as follows: In Sec. II we ex-plain the theoretical foundations of SCDFT, includingSF and SOI, and in Sec. III, details of the mathematicalformulation and the implementation are shown. Next,we list the results together with the numerical conditionin Sec. IV, and present the study discussion in Sec. V.Finally, we summarize the study in Sec. VI. II. THEORY
In this section, we will explain in detail the SCDFTformulation including plasmon-aided mechanism , SFeffect , and SOI . We use the Hartree atomic unitsthroughout the paper. In this study, we only considerthe singlet superconductivity, while in Ref. 16, boththe singlet- and triplet-superconducting states were con-sidered. Within SCDFT, T c is obtained as a temper-ature where the following Kohn-Sham superconductinggap ∆ n k becomes zero at all the band n and wavenum-ber k :∆ n k = − (cid:88) n (cid:48) k (cid:48) K n k n (cid:48) k (cid:48) ( ξ n k , ξ n (cid:48) k (cid:48) )1 + Z n k ( ξ n k ) × ∆ n (cid:48) k (cid:48) (cid:112) ξ n (cid:48) k (cid:48) + ∆ n (cid:48) k (cid:48) tanh (cid:32) (cid:112) ξ n (cid:48) k (cid:48) + ∆ n (cid:48) k (cid:48) T (cid:33) , (1) a r X i v : . [ c ond - m a t . s up r- c on ] M a r where ξ n k is the Kohn-Sham eigenvalue measured fromthe Fermi level ( ε F ) at the band index n and wave-number k . ξ n k is obtained by solving the following spinorKohn-Sham equation: (cid:32) − ∇ + V KS ↑↑ ( r ) − ε F V KS ↑↓ ( r ) V KS ↓↑ ( r ) − ∇ + V KS ↓↓ ( r ) − ε F (cid:33) (cid:18) ϕ n k ↑ ( r ) ϕ n k ↓ ( r ) (cid:19) = ξ n k (cid:18) ϕ n k ↑ ( r ) ϕ n k ↓ ( r ) (cid:19) , (2)where ϕ n k σ ( r ) is the σ component of the spinor Kohn-Sham orbital at ( n, k ), and V KS σσ (cid:48) ( r ) is the σσ (cid:48) compo-nent of the Kohn-Sham potential with SOI ( σ, σ (cid:48) = ↑ , ↓ ).Due to the off-diagonal part of the Kohn-Sham potential,the spin-up state and the spin-down state are hybridized.Therefore, the Kohn-Sham eigenvalue ξ n k does not havea spin index ( σ ). The non-linear gap equation (1) shouldbe solved numerically at each temperature. The integra-tion kernel K n k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) indicates the superconducting-pair breaking and creating interaction and comprises thefollowing three terms: K n k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) ≡ K epn k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) + K een k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) )+ K sfn k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) , (3)namely, the electron-phonon, Coulomb repulsion, andspin-fluctuation kernel, respectively. However, the renor-malization factor Z n k ( ξ n k ) comprises only the electron-phonon and spin-fluctuation terms as follows. Z n k ( ξ ) ≡ Z epn k ( ξ ) + Z sfn k ( ξ ) , (4)because the Coulomb-repulsion contribution to this fac-tor is already included in the Kohn-Sham eigenvalue ξ n k .The temperature T is defined by considering the Boltz-mann constant k B = 1.Let us explain each term in the kernel and the renor-malization factor below. The electron-phonon kernel K ep and renormalization factor Z ep are given by K epn k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) = 2tanh[ ξ/ (2 T )] tanh[ ξ (cid:48) / (2 T )] (cid:88) ν | g νn k n (cid:48) k (cid:48) | × [ I ( ξ, ξ (cid:48) , ω k (cid:48) − k ν ) − I ( ξ, − ξ (cid:48) , ω k (cid:48) − k ν )] , (5) Z epn k ( ξ ) = − ξ/ (2 T )] (cid:88) n (cid:48) k (cid:48) ν | g νn k n (cid:48) k (cid:48) | × [ J ( ξ, ξ n (cid:48) k (cid:48) , ω k (cid:48) − k ν ) + J ( ξ, − ξ n (cid:48) k (cid:48) , ω k (cid:48) − k ν )] , (6)where ω q ν is the phonon frequency at wave-number q and branch ν . I ( ξ, ξ (cid:48) , ω ) and J ( ξ, ξ (cid:48) , ω ) are derived withthe Kohn-Sham perturbation theory , and are writtenas follows : I ( ξ, ξ (cid:48) , ω ) = f T ( ξ ) f T ( ξ (cid:48) ) n T ( ω ) × (cid:34) e ξ/T − e ( ξ (cid:48) + ω ) /T ξ − ξ (cid:48) − ω − e ξ (cid:48) /T − e ( ξ + ω ) /T ξ − ξ (cid:48) + ω (cid:35) , (7) J ( ξ, ξ (cid:48) , ω ) = ˜ J ( ξ, ξ (cid:48) , ω ) − ˜ J ( ξ, ξ (cid:48) , ω ) , (8) ˜ J ( ξ, ξ (cid:48) , ω ) = − f T ( ξ ) + n T ( ξ ) ξ − ξ (cid:48) − ω × (cid:20) f T ( ξ (cid:48) ) − f T ( ξ − ω ) ξ − ξ (cid:48) − ω − f T ( ξ − ω ) f T ( − ξ (cid:48) + ω ) T (cid:21) , (9)where f T ( ξ ) and n T ( ω ) are the Fermi-Dirac and theBose-Einstein distribution function, respectively. Thefunctions I ( ξ, ξ (cid:48) , ω ) and J ( ξ, ξ (cid:48) , ω ) yield a temperature-dependent retardation effect. The electron-phonon ver-tex g between Kohn-Sham orbitals indexed with ( n, k )and ( n (cid:48) , k + q ), and the phonon ( q , ν ) is computed as g νn k n (cid:48) k + q = (cid:90) d r (cid:88) σσ (cid:48) = ↑ , ↓ ϕ ∗ n (cid:48) k + q σ ( r ) ϕ n k σ (cid:48) ( r ) × (cid:88) τ η τ q ν · δ q τ V KS σσ (cid:48) ( r ) (cid:112) M τ ω q ν , (10)where M τ is the mass of atom labeled by τ , η τ q ν is thepolarization vector of phonon ( q , ν ) and atom τ , and δ q τ V KS σσ (cid:48) ( r ) is the Kohn-Sham potential deformed by theperiodic displacement of atom τ and wave number q δ q τ V KS σσ (cid:48) ( r ) = (cid:88) R e i q · R δV KS σσ (cid:48) [ { r τ R } ]( r ) δ r τ R , (11)where r τ R is the position of the atom τ at the cell R . Wehave obtained the deformation potential by the phononcalculation, based on density functional perturbation the-ory (DFPT) . The electron-phonon kernel K ep is alwaysnegative; therefore, it makes a positive contribution informing the Cooper pair. However, the electron-phononrenormalization factor Z ep weakens the effect caused bythe kernels.The electron-electron repulsion kernel K ee in Eq. (3)is K een k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) = 2 π (cid:90) ∞ dω | ξ | + | ξ (cid:48) | ( | ξ | + | ξ (cid:48) | ) + ω V een k n (cid:48) k (cid:48) ( iω ) , (12)where V een k n (cid:48) k (cid:48) ( iω ) is the dynamically screened exchangeintegral between the Kohn-Sham orbitals ( n, k ) and( n (cid:48) , k (cid:48) ) V een k n (cid:48) k (cid:48) ( iω ) = (cid:90) (cid:90) d rd r (cid:48) V RP A ( r , r (cid:48) , iω ) × ρ (0) n k n (cid:48) k (cid:48) ( r ) ρ (0) ∗ n k n (cid:48) k (cid:48) ( r (cid:48) ) , (13) ρ (0) n k n (cid:48) k (cid:48) ( r ) = (cid:88) σ = ↑ , ↓ ϕ ∗ n k σ ( r ) ϕ n (cid:48) k (cid:48) σ ( r ) . (14)In this study, we have computed the screened Coulombinteraction V RP A by applying the random phase approx-imation (RPA) as V RP A ( r , r (cid:48) , iω ) = 1 | r − r (cid:48) | + (cid:90) (cid:90) d r d r × V RP A ( r , r , iω )Π ( r , r , iω ) 1 | r − r (cid:48) | , (15)where Π is the electronic susceptibility of the Kohn-Sham system (the non-perturbed susceptibility). Thiselectronic susceptibility is the α = 0 part of the followingsusceptibilities of the Kohn-Sham system:Π αα KS ( r , r (cid:48) , iω ) = (cid:88) kk (cid:48) nn (cid:48) θ ( − ξ n k ) − θ ( − ξ n (cid:48) k (cid:48) ) ξ n k − ξ n (cid:48) k (cid:48) + iω × ρ ( α ) n k n (cid:48) k (cid:48) ( r ) ρ ( α ) ∗ n k n (cid:48) k (cid:48) ( r (cid:48) ) , (16)where α takes 0, x , y , z , and ρ ( x ) n k n (cid:48) k (cid:48) ( r ) = (cid:88) σ = ↑ , ↓ ϕ ∗ n k σ ( r ) ϕ n (cid:48) k (cid:48) − σ ( r ) , (17) ρ ( y ) n k n (cid:48) k (cid:48) ( r ) = (cid:88) σ = ↑ , ↓ σϕ ∗ n k σ ( r ) ϕ n (cid:48) k (cid:48) − σ ( r ) , (18) ρ ( z ) n k n (cid:48) k (cid:48) ( r ) = (cid:88) σ = ↑ , ↓ σϕ ∗ n k σ ( r ) ϕ n (cid:48) k (cid:48) σ ( r ) . (19)The spin susceptibility Π xx KS , Π yy KS , and Π zz KS are used inthe spin-fluctuation term later on. Because of the factor[ θ ( − ξ n k ) − θ ( − ξ n (cid:48) k (cid:48) )] / ( ξ n k − ξ n (cid:48) k (cid:48) + iω ), these suscepti-bilities are affected largely by the electronic states in thevicinity of Fermi surfaces.We propose the spin-fluctuation (SF) kernel K sf inEq. (3) and the renormalization Z sf in Eq. (4) con-structed using the noncollinear spinor wavefunctions.The following formulation is an extension of those quan-tities in the collinear magnetism . K sfn k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) = 2 π (cid:90) ∞ dω | ξ | + | ξ (cid:48) | ( | ξ | + | ξ (cid:48) | ) + ω × Λ sfn k n (cid:48) k (cid:48) ( iω ) , (20) Z sfn k ( ξ ) = 1 π (cid:88) n (cid:48) k (cid:48) (cid:90) ∞ dω ( | ξ | + | ξ n (cid:48) k (cid:48) | ) − ω [( | ξ | + | ξ n (cid:48) k (cid:48) | ) + ω ] × Λ sfn k n (cid:48) k (cid:48) ( iω ) , (21)whereΛ sfn k n (cid:48) k (cid:48) ( iω ) = (cid:88) α = x,y,z (cid:90) (cid:90) d rd r (cid:48) Λ sfαα ( r , r (cid:48) , iω ) × ρ ( α ) n k n (cid:48) k (cid:48) ( r ) ρ ( α ) ∗ n k n (cid:48) k (cid:48) ( r (cid:48) ) . (22)Λ sfn k n (cid:48) k (cid:48) has a similar form to the screened exchange in-tegral of Eq. (13), and it involves the summation overthe x, y , and z components of the following SF-mediatedinteraction:Λ sfαα ( r , r (cid:48) , iω ) = − (cid:90) (cid:90) d r d r × I ααXC ( r , r )Π αα ( r , r , iω ) I ααXC ( r , r (cid:48) ) , (23)where Π αα is the spin susceptibility of the interactingsystem as Π αα ( r , r (cid:48) , iω ) = Π αα KS ( r , r (cid:48) , iω ) + (cid:90) (cid:90) d r d r × Π αα ( r , r , iω ) I ααXC ( r , r )Π αα KS ( r , r (cid:48) , iω ) . (24)In Eqs. (23) and (24), the spin-spin interaction is in-cluded through the exchange correlation kernel: I ααXC ( r , r (cid:48) ) ≡ δ E XC δm α ( r ) δm α ( r (cid:48) ) (25)which is the second-order functional derivative of theexchange-correlation energy E XC with respect to the spindensity along the α direction, m α . We have used theresults of the standard density functional calculationsof the normal (non-superconducting) state to calculatethe abovementioned quantities. Therefore, we have com-puted T c by solving the gap equation (1) as a post-processof the calculations of the normal state. The treatmentdescribed, known as the decoupling approximation, isknown to be reliable when the bandwidth and the su-perconducting gap energy scales are largely different . III. IMPLEMENTATION
In this section, we will explain the practical procedureto perform the calculations explained in the previous sec-tion.
A. Evaluation of exchange integrals with Fouriertransformation
The exchange integrals with Coulomb interaction part V een k n (cid:48) k (cid:48) of Eq. (13) and the SF part Λ sfn k n (cid:48) k (cid:48) of Eq. (22)can be computed efficiently using the Fourier transforma-tion as follows: First ρ ( α ) ( r ) in Eqs. (14) and (17)-(19)has the periodicity of the lattice vector R together withthe phase factor from the Bloch theorem as ρ ( α ) n k n (cid:48) k (cid:48) ( r + R ) = e i ( k (cid:48) − k ) · R ρ ( α ) n k n (cid:48) k (cid:48) ( r ) . (26)Therefore ρ ( α ) ( r ) can be expanded with the Fourier com-ponents of the reciprocal lattice vectors G as ρ ( α ) n k n (cid:48) k + q ( r ) = (cid:88) G e i ( q + G ) · r ˜ ρ ( α ) n k n (cid:48) k + q ( G ) , (27)where ˜ ρ ( α ) ( G ) is defined by the Fourier transformationof ρ ( α ) ( r ) as follows:˜ ρ ( α ) n k n (cid:48) k + q ( G ) ≡ v uc (cid:90) uc d re − i ( q + G ) · r ρ ( α ) n k n (cid:48) k + q ( r ) . (28)Subsequently, the exchange integral of Eq. (13) is rewrit-ten as V een k n (cid:48) k + q ( iω ) = (cid:88) GG (cid:48) V q RP A ( G , G (cid:48) , iω ) × ˜ ρ (0) n k n (cid:48) k + q ( G )˜ ρ (0) ∗ n k n (cid:48) k + q ( G (cid:48) ) , (29) Structure/Charge-density optimizationwith medium-density k gridPhonon calculation (DFPT) with medium-density k gridLoop:coarse q gridCoulomb/SF interaction Charge density Electron-phonon coupling with coarse k gridOne-shot DFT calculation with dense k gridOne-shot DFT calculation with coarse k gridEnd loop:coarse q gridLoop:coarse q gridEnd loop:coarse q gridGap equation (SCDFT) with coarse k grid Solve gap equation (SCDFT) at T =0Solve gap equation at T Gap averagedover Fermi surfaces hj ¢ ji < ¡ [meV] ? T c min = 0 T c max = T Solve gap equation at T =( T c min + T c max )/2Averaged gap hj ¢ ji < ¡ ¢ ? ¢ ´hj ¢ ji T =2 ¢ /3.54 YesAveraged gap hj ¢ ji < ¡ ¢ ?Loop:10 timesEnd loop:10 times T c min = T T c max = TT c = 0Yes NoYesNo T ( TT c =( T c min + T c max )/2(1)(2a)(3a)(3b)(3c)(4)(2b) In this step,This quantity is used that quantity is computed.to compute that quantity.OutputLegend No (a) (b) Eq. (11)Eq. (10) Eq. (32)Eq. (29)Eq. (36) Eq. (35)
FIG. 1. (a) The flow chart to perform the SCDFT calculation for each material in this study. The indices of steps are thesame as those mentioned in Sec. III D. (b) The flow chart of the bisection method used to compute T c . where V q RP A is the Fourier component of the screenedCoulomb interaction as follows: V q RP A ( G , G (cid:48) , iω ) ≡ (cid:90) (cid:90) d rd r (cid:48) e i ( q + G ) · r e − i ( q + G (cid:48) ) · r (cid:48) V RP A ( r , r (cid:48) , iω ) . (30)However, the values of V RP A ( r , r , iω ) need not be found,since we show that V q RP A ( G , G (cid:48) , iω ) at each q can becomputed separately using the Bloch theorem as shownbelow. By substituting V RP A ( r , r (cid:48) , iω ) of Eq. (15) intoEq. (30), we obtain V q RP A ( G , G (cid:48) , iω ) = 4 πδ GG (cid:48) | q + G | + (cid:88) G V q RP A ( G , G , iω )Π q KS ( G , G (cid:48) , iω ) 4 π | q + G (cid:48) | = (cid:20) | q + G (cid:48) | δ GG (cid:48) π − Π q KS ( G , G (cid:48) , iω ) (cid:21) − , (31)where Π αα q KS ( G , G (cid:48) , iω ) is the Fourier component of thesusceptibilities of the Kohn-Sham system of Eq. (16)given byΠ αα q KS ( G , G (cid:48) , iω ) ≡ (cid:90) (cid:90) d rd r (cid:48) e i ( q + G ) · r e − i ( q + G (cid:48) ) · r (cid:48) Π αα KS ( r , r (cid:48) , iω )= (cid:88) k nn (cid:48) θ ( − ξ n k ) − θ ( − ξ n (cid:48) k + q ) ξ n k − ξ n (cid:48) k + q + iω ˜ ρ ( α ) n k n (cid:48) k + q ( G )˜ ρ ( α ) ∗ n k n (cid:48) k + q ( G (cid:48) ) . (32)In the derivation of Eq. (31), we used the periodicity ofeach term with respect to the lattice vector. The factor[ θ ( − ξ n k ) − θ ( − ξ n (cid:48) k + q )] / ( ξ n k − ξ n (cid:48) k + q + iω ) in the sus-ceptibilities varies rapidly in the vicinity of Fermi sur-faces, and we need a dense k grid to compute it accu-rately, which may require a huge numerical cost. There-fore, we use the reverse interpolation scheme explained inSec. III.C.1 of Ref. 23. In this scheme, we compute the ex-plicitly energy-dependent factor with a dense k grid whilewe compute the other parts using a coarse k grid because˜ ρ ( α ) n k n (cid:48) k + q ( G )˜ ρ ( α ) ∗ n k n (cid:48) k + q ( G (cid:48) ) varies more smoothly than theenergy-dependent factor. The SF term can be computedin the same manner at each q separately. The Fouriercomponent of the SF-mediated interaction of Eq. (23) isΛ sf,αα q ( G , G (cid:48) , iω ) = − (cid:88) G , G I αα q XC ( G , G ) × (cid:104) (Π αα q KS ( G , G )) − − I αα q XC ( G , G ) (cid:105) − I αα q XC ( G , G (cid:48) ) , (33)where I αα q XC is the Fourier component of the exchange cor-relation kernel of Eq. (25). In this study, we employed thelocal density approximation (LDA) to describe this ker-nel as follows: We approximate the exchange-correlationenergy as E LDAXC = (cid:90) d rε homXC ( ρ ( r ) , | m ( r ) | ) ρ ( r ) , (34)where ρ ( r ) is the electronic charge density and ε homXC ( ρ, m )is the exchange-correlation energy density of the homo-geneous electon gus whose charge and spin density are ρ and m . The exchange-correlation kernel in Eq. (25)becomes I LDA,ααXC ( r , r (cid:48) )= δ ( r − r (cid:48) ) ρ ( r ) ∂ ε ( ρ ( r ) , | m | ) ∂ | m | ∂ | m | (cid:18) m α | m | (cid:19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m = m ( r ) . (35)Since we perform the non-magnetic calculation in thisstudy, we take the m ( r ) → I αα q XC ( G , G (cid:48) ) does not depend on q . Thisis equivalent to the adiabatic local density approxi-mation (ALDA) in time-dependent density functionaltheory . B. Auxiliary gap equation
We have solved the gap equation (1) using the auxiliaryenergy axis to capture the rapid change of the explic-itly energy-dependent function in the vicinity of Fermi surfaces. In this method, the gap function ∆ n k dependsalso on the auxiliary energy; the auxiliary gap function˜∆ n k ( ξ ) satisfies ˜∆ n k ( ξ n k ) = ∆ n k . Subsequently, the gapequation (1) becomes˜∆ n k ( ξ ) = − (cid:90) ∞−∞ dξ (cid:48) (cid:88) n (cid:48) k (cid:48) D n (cid:48) k (cid:48) ( ξ (cid:48) ) K n k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) )1 + Z n k ( ξ ) × ˜∆ n (cid:48) k (cid:48) ( ξ (cid:48) ) (cid:113) ξ (cid:48) + ˜∆ n (cid:48) k (cid:48) ( ξ (cid:48) ) tanh (cid:113) ξ (cid:48) + ˜∆ n (cid:48) k (cid:48) ( ξ (cid:48) )2 T , (36)where D n k ( ξ ) is the ( n, k )-resolved density of states. Inthe same manner, the electron-phonon and SF renormal-ization factors of Eqs. (6) and (21) become Z epn k ( ξ ) = − ξ/ (2 T )] (cid:90) ∞−∞ dξ (cid:48) (cid:88) n (cid:48) k (cid:48) ν D n (cid:48) k (cid:48) ( ξ (cid:48) ) | g νn k n (cid:48) k (cid:48) | × [ J ( ξ, ξ (cid:48) , ω k (cid:48) − k ν ) + J ( ξ, − ξ (cid:48) , ω k (cid:48) − k ν )] (37)and Z sfn k ( ξ ) = 1 π (cid:90) ∞−∞ dξ (cid:48) (cid:88) n (cid:48) k (cid:48) D n (cid:48) k (cid:48) ( ξ (cid:48) ) (cid:90) ∞ dω × ( | ξ | + | ξ n (cid:48) k (cid:48) | ) − ω [( | ξ | + | ξ n (cid:48) k (cid:48) | ) + ω ] Λ sfn k n (cid:48) k (cid:48) ( iω ) , (38)respectively, where we employ the reverse interpola-tion method again; the ( n, k )-resolved density of states D n k ( ξ ) is computed with the dense k grid, while theother parts are computed on the coarse k grid; finally,we combine the parts yielded. C. Frequency integral
The integration in Eq. (12) involves the frequency ω spanning [0, ∞ ]. Therefore, to perform the integrationnumerically, we change the variable as follows: ω = ( | ξ | + | ξ (cid:48) | ) 1 + x − x . (39)Then Eq. (12) becomes K een k n (cid:48) k (cid:48) ( ξ, ξ (cid:48) ) = 2 π (cid:90) − dx
11 + x × V een k n (cid:48) k (cid:48) (cid:18) i ( | ξ | + | ξ (cid:48) | ) 1 + x − x (cid:19) . (40)When | ξ | + | ξ (cid:48) | = 0, this integration becomes V een k n (cid:48) k (cid:48) (0).To obtain V een k n (cid:48) k (cid:48) ( iω ) at an arbitrary ω , we first com-pute V een k n (cid:48) k (cid:48) ( iω ) at discrete non-uniform points ω =0 , ω , ω , · · · , ∞ and interpolate them. Using the sametransformation as in Eq. (39), the frequency integral inthe SF renormalization term (38) is performed as follows: Z sfn k ( ξ ) = 1 π (cid:90) ∞−∞ dξ (cid:48) (cid:88) n (cid:48) k (cid:48) D n (cid:48) k (cid:48) ( ξ (cid:48) ) (cid:90) − dx − x ( | ξ | + | ξ n (cid:48) k (cid:48) | )(1 + x ) × Λ sfn k n (cid:48) k (cid:48) (cid:18) i ( | ξ | + | ξ (cid:48) | ) 1 + x − x (cid:19) . (41)The numerical integration with respect to the variable x ranging [-1,1] can be performed using the Gauss quadra-ture. D. Overall procedure
Figure 1 (a) shows the calculation flow. We employ thefollowing three different wavenumber grids to efficientlyperform the Brillouin-zone integrals: coarse grid:
To reduce the computational cost, we usea coarse grid for the wavenumber q of phonons andsusceptibilities. The grid is shifted by half of itsstep to avoid the singularity at the Γ point. Thisgrid is also used for solving the gap equation. medium-density grid: The atomic structure and thecharge density are optimized with the self-consistent field calculation using k grid denser thanthe coarse grid. This k grid is used to prepare elec-tronic states in the DFPT calculation. dense grid: To treat the explicitly energy-dependentfactor in the calculations of susceptibilities and the( n, k )-dependent density of states in the gap equa-tion (36), a dense k grid is employed.The overall calculations are performed as follows:1. First, we optimize the atomic structure and thecharge density using the standard density func-tional calculation with the medium-density k grid.The following calculation is performed on this opti-mized atomic structure and with the charge density.2. We compute the electron-phonon interaction andthe frequency of phonons, whose wavenumber q ison the half-grid shifted coarse grid. This step isfurther split into the following two sub-steps:(a) The phonon calculation based on DFPT isperformed; the electronic states used in thiscalculation have a wavenumber k on themedium-density grid.(b) Electron-phonon vertex of Eq. (10) betweenthe Kohn-Sham orbitals ( n, k ) and ( n (cid:48) , k + q )is computed, where the wavenumber k is onthe coarse grid.3. Next, we compute the exchange integrals ofscreened Coulomb and SF-mediated interactionwhose transitional momentum q is on the half-gridshifted coarse grid. This step is split into the fol-lowing three sub-steps: (a) One-shot DFT calculation on the dense k gridis performed. The resulting energy disper-sion ξ n k is later used to compute the explicitlyenergy-dependent term in the susceptibilitiesin Eq. (32).(b) One-shot DFT calculation on the coarse k grid with and without a half-grid shift is per-formed. These two grids are connected by thetransitional momentum q on the coarse a half-grid shifted grid. The resulting Kohn-Shamorbitals will be used to compute ρ ( α ) ( r ) of Eqs.(14) and (17)-(19).(c) The exchange integrals of screened Coulomband SF-mediated interaction between theKohn-Sham orbitals ( n, k ) and ( n (cid:48) , k + q ) arecomputed, where the wavenumber k is on thecoarse grid.4. The gap equation within SCDFT is solved on thecoarse k grid at each temperature. Then T c is ob-tained as a minimum temperature where all ∆ n k ( ξ )vanish.The transition temperature T c is found using the bisec-tion method explained in Fig. 1 (b). While the initiallower limit of T c is set to zero, the initial upper limit isset to the T c estimated by the Bardeen-Cooper-Schrieffertheory (2∆ / .
54, where ∆ is the superconductinggap averaged over Fermi surfaces at zero kelvin). If thereis a finite gap even at this upper limit, although it rarelyoccurs, we double the initial upper limit; then we repeatthe bisection step ten times and find T c . IV. RESULT
In this section we will first explain the numerical condi-tion of this study, then show the result of the benchmark.The numerical condition is as follows: We use the DFTcode
Quantum ESPRESSO which employs planewaves and pseudopotentials. Perdew-Burke-Ernzerhof’sdensity functional based on generalized gradient ap-proximation (GGA) is used. We use the optimizednorm-conserving pseudopotential library provided bySchlipf-Gygi (SG15) . The energy cutoff for the wavefunctions of each element is specified using the criteriaand the convergence profiles in the Standard Solid StatePseudopotentials . We are using the optimized tetra-hedron method to perform the Brillouin-zone integra-tion. The number of grids along each reciprocal latticevector is proportional to the length of that vector. Ta-ble I presents the above explained conditions for eachelement. The SCDFT calculation is performed using the Superconducting Toolkit software package. Wehave set the medium-density grid twice the density ofthe coarse grid and set the dense grid twice the densityof the medium-density grid; for example, 8 , 16 , and32 grids were used for the coarse, medium-density, and TABLE I. Structure, cutoff for the plane wave for wave func-tions, q -grid for phonon, and the error of the lattice constant (cid:104) ∆ a/a exp (cid:105) ≡ ( V calc /V exp ) / − V calc and V exp are the cal-culated and experimental unit-cell volume, respectively.structure cutoff [Ry] coarse grid (cid:104) ∆ a/a exp (cid:105) [%]Be hcp 65 10 × × × × × × × × × × × × × × × × × × × × × × α -Ga 150 5 × × × × × × × × × × × × × × × × × × × × × × × × × × × × β -Sn 50 7 × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × dense grids, respectively; all for Al. The minimum scaleand the number of points of the non-uniform auxiliaryenergy grid used in the gap equation (36) were set to10 − Ry and 100, respectively. In the calculation of themagnetic exchange-correlation kernel of Eq. (25), we ig-nore the gradient correction. For stabilizing the phononcalculation, the lattice constants and the internal atomiccoordinates are optimized; the deviation between the op-timized and the experimental lattice constants are listedin Tbl. I. To compute the susceptibilities in Eq. (32) andsolve the gap equation (36), we have included 40 × N atom (20 × N atom ) empty bands for the calculation with (with-out) SOI, where N atom is the number of atoms per unit cell.Next, we move onto the results. Figure 2 shows theexperimental T c ( T exp c ) , the theoretical T c computedwith and without SOI/SF in a periodic-table form. Also,to examine the effects of the electron-phonon interac-tion, the screened Coulomb repulsion, and the spin-fluctuation, we are showing the following quantities in thesame figure: The density of states (DOS) at the Fermilevel divided by the number of atoms affects the strengthof the mean-field. The Fr¨ohlich’s mass-enhancement pa-rameter λ = (cid:88) q ν λ q ν , (42)and the averaged phonon frequencies ω ln = exp (cid:34) λ (cid:88) q ν λ q ν ln( ω q ν ) (cid:35) , (43)appear in the conventional McMillan formula T c = ω ln . (cid:20) − . λ ) λ − µ ∗ (1 + 0 . λ ) (cid:21) (44)which has been used to estimate T c semi-empirically withan adjustable parameter µ ∗ . The ( q , ν )-dependent mass-enhancement parameter λ q ν is computed as follows: λ q ν = 2 D (0) ω q ν (cid:88) k nn (cid:48) | g νn k n (cid:48) k + q | δ ( ξ n k ) δ ( ξ n (cid:48) k + q ) , (45)where D (0) is the density of states at the Fermi level.We are calculating here the Brillouin-zone integral, in-cluding two delta functions, using the dense k grid to-gether with the optimized tetrahedron method . Be-cause of the double delta function δ ( ξ n k ) δ ( ξ n (cid:48) k + q ), thissummation involves the electron-phonon vertices only be-tween the electronic states at the Fermi level. There-fore, the Fr¨ohlich’s mass-enhancement parameter λ inEq. (42) indicates the retarded phonon-mediated inter-action (2 | g | /ω ) averaged over Fermi surfaces times thedensity of state at the Fermi level. Similarly, ω ln indicatesthe typical frequency of phonons which couples largelywith the electronic states at the Fermi level. Therefore, λ and ω ln closely relate to the electron-phonon contribu-tion to T c . In an analogous fashion to λ , we are showingparameters for the Coulomb repulsion and SF as µ C = 1 D (0) (cid:88) kk (cid:48) nn (cid:48) K een k n (cid:48) k (cid:48) δ ( ξ n k ) δ ( ξ n (cid:48) k (cid:48) ) , (46)and µ s = 1 D (0) (cid:88) kk (cid:48) nn (cid:48) K sfn k n (cid:48) k (cid:48) δ ( ξ n k ) δ ( ξ n (cid:48) k (cid:48) ) , (47)respectively. These parameters are the Coulomb[Eq.(12)] and SF [Eq.(20)] kernels averaged over Fermisurfaces times the density of states at the Fermi level. . . .
527 67 1 .
04 0 .
224 0 . .
032 0 .
013 0 .
062 685 0 .
161 0 .
080 0 .
017 0 0 0 .
472 139 0 .
181 0 .
332 0 . .
428 253 0 .
237 0 .
280 0 . . .
89 0 .
404 302 0 .
402 0 .
251 0 .
076 0 0 0 .
768 72 0 .
132 0 .
347 0 . .
58 148 0 .
131 0 .
484 0 . . .
01 205 0 .
398 0 .
523 0 . . .
593 0 .
894 207 0 .
498 0 .
258 0 . . .
77 231 1 .
26 0 .
512 0 . C r : M n : F e : C o : N i : .
895 39 0 .
159 0 .
357 0 . . .
51 0 .
827 270 0 .
472 0 .
278 0 . . .
35 240 0 .
386 0 .
461 0 . .
275 138 0 .
137 0 .
162 0 . . . .
236 48 0 .
702 0 .
143 0 . . .
929 193 0 .
967 0 .
286 0 . . . .
207 112 0 . .
152 0 . N o d a t a du e t o i m ag i n a r y ph o n o n f r e q u e n c i e s . . .
30 149 0 .
923 0 .
396 0 . .
26 0 .
18 0 .
432 225 0 .
265 0 .
151 0 . . . .
777 177 0 .
635 0 .
253 0 . . .
61 0 .
606 228 0 .
409 0 .
216 0 . .
73 145 0 .
333 0 .
977 0 . . . .
497 86 0 .
867 0 .
234 0 .
051 0 .
13 0 1 .
19 63 0 .
205 0 .
286 0 . . . .
445 47 1 .
03 0 .
215 0 . .
051 0 .
292 102 0 .
186 0 .
129 0 . .
26 112 0 .
528 0 .
798 0 . .
83 0 .
26 0 .
928 203 0 .
422 0 .
339 0 . .
38 24 0 .
155 0 .
429 0 . . .
48 154 1 .
24 0 .
429 0 . . . .
593 265 0 .
438 0 .
190 0 . .
10 0 2 .
15 101 0 .
068 0 .
622 0 . . . .
467 58 1 .
17 0 .
229 0 . . . .
234 120 0 .
573 0 .
157 0 . .
017 0 0 .
305 216 0 .
119 0 .
163 0 . . . .
207 112 0 .
509 0 .
153 0 .
027 3 . . .
235 49 0 .
681 0 .
142 0 .
025 2 . . .
402 47 0 .
752 0 .
197 0 . . . .
564 58 1 .
45 0 .
241 0 . . . .
508 87 0 .
956 0 .
238 0 . . . .
233 119 .
575 0 .
159 0 . . .
88 0 .
404 303 0 .
401 0 .
251 0 .
076 0 0 0 .
276 139 0 .
137 0 .
161 0 .
022 0 .
10 0 .
081 0 .
309 98 0 .
199 0 .
136 0 . .
306 216 0 .
119 0 .
163 0 .
021 1 . .
26 1 .
84 125 0 .
491 0 .
668 0 . .
74 148 0 .
325 0 . - . .
95 0 1 .
34 240 0 .
384 0 .
458 0 .
269 0 .
85 0 .
26 0 .
941 202 0 .
436 0 .
342 0 . . .
50 0 .
613 230 0 .
395 0 .
221 0 . . . .
662 189 0 .
482 0 .
219 0 . . .
49 0 .
831 270 0 .
479 0 .
280 0 . . .
910 201 0 .
890 0 .
280 0 . . . .
590 265 0 .
443 0 .
190 0 .
057 0 .
094 0 .
058 0 .
381 227 0 .
237 0 .
137 0 . . .
48 155 1 .
23 0 .
431 0 .
203 6 . . .
34 149 0 .
912 0 .
408 0 . .
032 0 0 .
472 140 0 .
180 0 .
332 0 .
213 0 0 0 .
768 72 0 .
133 0 .
347 0 .
270 0 0 0 .
896 39 0 .
158 0 .
357 0 .
280 0 0 1 .
38 23 0 .
159 0 .
430 0 . .
078 0 0 .
428 253 0 .
237 0 .
279 0 . No data due to lack of pseudopotential . .
59 0 .
894 207 0 .
498 0 .
258 0 . . .
77 231 1 .
27 0 .
512 0 .
722 0 0 1 .
96 100 0 .
093 0 .
567 0 . .
58 148 0 .
131 0 .
484 0 . . .
01 205 0 .
398 0 .
521 0 . L i : . B e : . A l : . N a : M g : C a : K : R b : S r : Y : Z r : . N b : . M o : . T c : . R u : . R h : . P d : A g : C d : . I n : . C s : B a : L a : . H f : . T a : . W : . R e : . O s : . I r : . P t : A u : H g : . T l : . P b : . S c : Z n : . T i : . V : . G a : . Sn : . C u E l e m e n t : T c ( E x p . ) [ K ] T c ( w o S F ) T c ( w S F ) D O S ( " F ) [ e V ¡ ] ! l n [ K ] ¸ ¹ C ¹ s T c ( w o S F ) T c ( w S F ) D O S ( " F ) [ e V ¡ ] ! l n [ K ] ¸ ¹ C ¹ s w o S O w S O No data due to imaginary phonon frequency S k i p b ec a u s e o f c o m p li c a t e d s t r u c t u r e M ag n e t i c o r d e r M ag n e t i c o r d e r M ag n e t i c o r d e r M ag n e t i c o r d e r M ag n e t i c o r d e r No data due to lack of pseudopotential N o d a t a du e t o i m ag i n a r y ph o n o n f r e q u e n c i e s N o d a t a du e t o i m ag i n a r y ph o n o n f r e q u e n c i e s N o d a t a du e t o i m ag i n a r y ph o n o n f r e q u e n c i e s N o d a t a du e t o i m ag i n a r y ph o n o n f r e q u e n c i e s F I G . . T h ee x p e r i m e n t a l T c ( T e x p c ) a nd t h e t h e o r e t i c a l T c c o m pu t e d w i t h a nd w i t h o u t S O I / S F . I n a dd i t i o n , w e s h o w t h e d e n s i t y o f s t a t e s ( D O S ) d i v i d e db y t h e nu m b e r o f a t o m s , t h e a v e r ag e dph o n o n f r e q u e n c i e s ω l n , t h e F r ¨o h li c h ’ s m a ss - e nh a n c e m e n t p a r a m e t e r λ , t h e a v e r ag e d C o u l o m b i n t e r a c t i o n µ C i n E q . ( ) , a nd t h e a v e r ag e dS F t e r m µ s i n E q . ( ) . C y a n , g r ee n , r e d , b l u e , a nd y e ll o w c e ll s i nd i c a t e t h e a l k a li n e ( e a r t h ) m e t a l s , t h e s p - o r b i t a l m e t a l s , t h e g r o up m e t a l s , t h e n o b l e m e t a l s , a nd t h e o t h e r t r a n s i t i o n m e t a l s , r e s p e c t i v e l y . B e N a 2 M g 13 A l K C a 3 S c T i V C u Z n G a 1 R b S r Y Z r N b M o 7 T c R u R h P d A g 12 C d I n Sn C s B a 3 L a 4 H f T a 6 W R e O s I r P t A u H g 13 T l P b ° S [ m J K - m o l - ] SOI :no SOI :Exp. :
FIG. 3. Theoretical and experimental Sommerfeld coefficient γ S . The horizontal axis is the atomic symbol together withthe group of the periodic table. Blue crosses and green circles indicate the Sommerfeld coefficient computed with and withoutSOI while magenta “+” with an error bar indicates the experimental value of the Sommerfeld coefficient. B e N a 2 M g 13 A l K C a 3 S c T i V C u Z n G a 1 R b S r Y Z r N b M o 7 T c R u R h P d A g 12 C d I n Sn C s B a 3 L a 4 H f T a 6 W R e O s I r P t A u H g 13 T l P b T c [ K ] SOI+SF :SOI :SF :no SOI/SF :Exp. : 20 22
FIG. 4. Theoretical and experimental T c s. The horizontal axis is the atomic symbol together with the group of the periodictable. Downward (upward) triangles indicate the T c s computed with (without) SF. Filled (empty) triangles indicate the T c scomputed with (without) SOI. “+” indicates the experimental value of T c . The plot which ranges from zero to two Kelvin inthe upper panel is magnified into the bottom panel. We note that in Fig. (2), some results are absent be-cause of the following reason: Li has a vast unit cell ata low temperature, and it is computationally demand-ing. Cr, Mn, Fe, Co, and Ni show a magnetic order atthe low temperature. Therefore, the formalism of thespin-fluctuation (23) used in this study breaks down; the matrix δ GG (cid:48) − (cid:88) G Π αα q KS ( G , G ) I αα q XC ( G , G (cid:48) ) (48)does not become positive definite (the Stoner’s criterion)for those materials. For Be and Ba, there is no pseudopo-tential together with SOI in the pseudopotential libraryused in this study. Since we are trying to unify the condi-0tion of the calculation for all elements, we leave the resultof these two materials together with SOI blank. For Y,Zr, In (with SOI), La, Hf, and Hg, we have obtainedimaginary phonon frequencies because of an artificiallong-range structure instability. For such cases, we couldhave not continue the calculation because of the break-down of the formulations of the electron-phonon kernel[Eq. (5)] and renormalization term [Eq. (6)]. Therefore,we leave the results for those cases blank.We have checked the convergences of the numericalresults with respect to the density of k , q , and the aux-iliary energy grids. For this purpose, we have selectedup the typical four materials, namely Al, V, Nb, andPb. Figure 5 shows the convergence of calculated T c (with SF, without SOI), the density of states at the Fermilevel, the averaged phonon frequencies ω ln , the Fr¨ohlich’smass-enhancement parameter λ , he averaged Coulombinteraction µ C in Eq. (46), and the averaged SF term µ s in Eq. (47) with respect to the density of the wave-number grids. Also, in the convergence check, we have setthe medium-density grid twice the density of the coarsegrid and have set the dense grid twice the density of themedium-density grid. Although the λ of V and Nb os-cillate within 4% and 7%, respectively, when we changethe density of the coarse grid, other quantities, including T c , are unchanged. We are showing the convergence of T c (with SF, without SOI) with respect to the numberof points for the energy integral in Eq. (36) in Fig. 6.This result shows the T c s have converged at the numeri-cal conditions described above.We compare λ in Eq. (42), ω ln in Eq. (43), and µ C inEq. (46) obtained in this study and earlier studies. Ta-ble II shows λ and ω ln in this study and earlier studies.Although there are small deviations due to the differ-ent exchange-correlation functional, these quantities areconsistent with those of earlier studies excepting the casefor Pb without SOI; the reported λ s of Pb without SOIare scattered and different from the experimental value( λ exp = 1 . . We also confirm that we reproducedthe averaged Coulomb interaction µ C in Eq. (46) in theearlier works for Al and Nb; µ C is 0.251 and 0.429 forAl and Nb, respectively in this study while those in theearlier studies are 0.236 and 0.488 .In Fig. 7, we have plotted the experimental T c togetherwith the T c calculated the most precisely in this work byincluding SOI and SF. We note that since there is nocomputed data for Be and In with SOI, we are showingthe data for them without SOI. Excepting Cd, Zn, and V,we have accurately reproduced experimental T c . Insidethe target elements of this benchmark, the three highest- T c elements observed experimentally are Nb (9.20 K),Tc (7.77 K), and Pb (7.19 K). Also, the three highest- T c elements in our calculation are Nb (7.470 K), Tc (6.019K), and Pb (6.010 K). Therefore, at least in the elementalmetals at ambient pressure, we can predict the highest- T c materials. T c S F [ K ] Al :V :Nb :Pb : 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 D O S [ e V ¡ /a t o m ]
50 100 150 200 250 300 ! l n [ K ] Al :V :Nb :Pb : 0.4 0.6 0.8 1 1.2 1.4 ¸ ¹ C Al :V :Nb :Pb : 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 6 ¹ s Coarse grid
FIG. 5. The convergence of calculated T c (with SF, withoutSOI), the density of states at the Fermi level, the averagedphonon frequencies ω ln , the Fr¨ohlich’s mass-enhancement pa-rameter λ , he averaged Coulomb interaction µ C in Eq. (46),and the averaged SF term µ s in Eq. (47) with respect to thedensity of the wave-number grids. We pick up Al, V, Nb, andPb. V. DISCUSSION
In this section, we discuss the results shown in the pre-vious section. To check the accuracy of the calculationof the electron-phonon interaction, we will plot in Fig. 3the calculated- and experimental- Sommerfeld coeffi-1
TABLE II. Fr¨ohlich’s mass-enhancement parameter λ in Eq. (42) and the averaged phonon frequencies ω ln in Eq. (43) computedin this study and earlier studies. Numerical conditions, the q grid of phonons, the k grid of electronic states for the calculationof phonon frequencies, the k grid and the smearing width for the integration in Eq. (45) are also written. In the Ref. 39 onlythe number of q and k points in the irreducible Brillouin-zone (IBZ) are provided. In this study and one of the earlier studies the tetrahedron method is employed for performing the Brillouin-zone integration in Eq. (45). λ ω ln [K] q grid k grid (phonon) k grid ( λ ) Smearing [eV] Functional Ref.Al (without SOI) 0.402 302 8 × × × ×
16 32 × ×
32 Tetrahedron GGA This work0.438 - 89 in IBZ - 1,300 in IBZ 0.272 LDA 390.417 314 8 × × × ×
12 32 × ×
32 0.272 LDA 150.44 270 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40V (without SOI) 1.26 231 9 × × × ×
18 36 × ×
36 Tetrahedron GGA This work1.19 245 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Cu (without SOI) 0.119 216 9 × × × ×
18 36 × ×
36 Tetrahedron GGA This work0.14 220 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Nb (without SOI) 1.24 154 8 × × × ×
16 32 × ×
32 Tetrahedron GGA This work1.26 185 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Mo (without SOI) 0.438 265 8 × × × ×
16 32 × ×
32 Tetrahedron GGA This work0.42 280 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Pd (without SOI) 0.333 145 8 × × × ×
16 32 × ×
32 Tetrahedron GGA This work0.35 180 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Ta (without SOI) 0.923 149 8 × × × ×
16 32 × ×
32 Tetrahedron GGA This work0.86 160 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 40Tl (without SOI) 1.03 47 6 × × × × × ×
12 Tetrahedron GGA This work1.0 - 8 × × × ×
16 36 × ×
24 0.2 LDA 19Tl (with SOI) 0.752 47 6 × × × × × ×
12 Tetrahedron GGA This work0.87 - 8 × × × ×
16 36 × ×
24 0.2 LDA 19Pb (without SOI) 1.04 67 6 × × × ×
12 24 × ×
24 Tetrahedron GGA This work1.20 - 89 in IBZ - 1,300 in IBZ 0.272 LDA 391.68 65 8 × × × ×
16 32 × ×
32 Tetrahedron LDA 401.08 - 8 × × × ×
16 32 × ×
32 0.2 LDA 191.24 - 8 × × × ×
16 40 × ×
40 0.272 LDA 3Pb (with SOI) 1.45 58 6 × × × ×
12 24 × ×
24 Tetrahedron GGA This work1.56 - 8 × × × ×
16 32 × ×
32 0.2 LDA 19 T c S F [ K ] Number of energy pointsAl :V :Nb :Pb :
FIG. 6. The convergence of the T c (with SF) with respectto the number of points for the energy integral in Eq. (36) forAl, V, Nb, and Pb. cient γ S which is the prefactor of the specific heat at a low temperature C v = γ S T + O ( T ) . (49)This coefficient is estimated using the density of statesand the Fr¨ohlich’s mass-enhancement parameter as fol-lows γ S = π D (0)3 (1 + λ ) . (50)The calculated γ S agrees very well with the experimen-tal value. For Re, Pt, and Pb, this agreement becomesimproved by including SOI; Since these heavy elementshave large SOI, this interaction is crucial to reproducethe experimental Sommerfeld coefficient.We have plotted the computed- and experimental- T c data contained in Fig. 2 into Fig. 4 to visualize the ef-fect by SOI and SF; we can detect the following trendsby inspecting this graph: SF always reduces T c s forthe elemental systems. This reduction becomes signifi-cant for the transition metals and is crucial to reproducethe experimental T c quantitatively. The mechanism ofthis reduction is explained as follows : For isotropicsuperconductors such as elemental materials, ferromag-netic spin-fluctuation becomes dominant and aligns the2 T c S O I + S F [ K ] T c exp [K]VCdZn Deviation = §
20 %
FIG. 7. Experimental and computed T c . For the theoreticalvalue, only the result computed with SOI and SF is repre-sented. The shaded region indicates that the deviation of T c ( T SOI+SF c /T exp c ) is less than 20 %. spin of electrons parallel. This effect breaks the singletCooper pair in the isotropic superconductors. On theother hand, in cuprates and iron-based superconductors,highly anisotropic antiferromagnetic spin-fluctuation be-comes dominant and enhances the superconducting gapswith sign-changes . In the transition metals, the effectof SF weakens with the increasing of the period numberin the periodic table. For example, µ s varies 0.722 (V) → → → → → → µ s of Pd becomes negative;this indicates the formulation of SF breaks down due tothe magnetic order. Although the magnetic order is anumerical artifact, this result shows that Pd has largerSF than that of Pt. This trend of SF can be explainedas follows : The electronic orbitals become delocalizedwith the increasing of the principal quantum number (3d → → µ s does not decrease with the increasingof the period number, i.e., this parameter varies 0.213(Na) → → → λ changes drastically by turning on the SOI. For Pb, -3-2-1 0 1 2 3 B e A l T i V Z n G a5 N b M o7 T c R u R h C d I n Sn T a6 W R e O s I r T l P b T c S O I + S F ¡ T c e x p [ K ] (a) T c S O I + S F / T c e x p (b) B e A l T i V Z n G a5 N b M o7 T c R u R h C d I n Sn T a6 W R e O s I r T l P b D O S [ e V ¡ /a t o m ] (c) ! l n [ K ] (d) B e A l T i V Z n G a5 N b M o7 T c R u R h C d I n Sn T a6 W R e O s I r T l P b ¸ (e) ¹ C (f) B e A l T i V Z n G a5 N b M o7 T c R u R h C d I n Sn T a6 W R e O s I r T l P b ¹ s (g) FIG. 8. We plot the following quantities for the elementalsystems which have finite T c : (a) The difference between theexperimental T c ( T exp c ) and theoretical T c computed with SOIand SF ( T SOI+SF c ). (b) The ratio between T exp c and T SOI+SF c .(c) The density of states at the Fermi level. (d) The averagedphonon frequencies ω ln . (e) The Fr¨ohlich’s mass-enhancementparameter λ . (f) The averaged Coulomb interaction µ C inEq. (46). (g) The averaged SF term µ s in Eq. (47). Thehorizontal axis is the atomic symbol together with the groupof the periodic table. We note that since there is no computeddata for Be and In with SOI, we show the data for themwithout SOI. this enhancement of λ (1.036 → ω ln decreases from 67 K to 58 K), the increased DOSat Fermi level (0.527 eV − → − ), and the en-hanced deformation potential δ q τ V KS ( r ) due to the SOIterm . These effects of SOI in Tc, Re, and Tl are op-posite to ones in Pb and Sn; SOI reduces the electron-phonon coupling as well as T c in these three materials.We can reproduce the absence of the superconductivityin alkaline, alkaline earth, and noble metals, exceptingPt and Au with SOI and SF; we have observed smallfinite T c for these two elements; we can reproduce thenon-superconductivity also in Sc by including SF whilewe observe T c = 2 .
711 K by ignoring SF. Since Sc hashighly localized 3d electrons, the SF largely reduces T c .For the group 12 elements (Zn and Cd), T c s are overesti-mated even if we include SF. For these materials, the SFeffect is small because the d orbitals are fully occupied.Finally, we have tried to find the factor which dom-inates the accuracy of T c . In Fig. 8, we have plottdthe following quantities for the elemental systems withfinite T c : (a) The difference between the experimental T c ( T exp c ) and theoretical T c computed with SOI and SF( T SOI+SF c ). (b) The ratio between T exp c and T SOI+SF c .(c) The density of states at the Fermi level. (d) Theaveraged phonon frequencies ω ln in Eq. (43). (e) TheFr¨ohlich’s mass-enhancement parameter λ in Eq. (42).(f) The averaged Coulomb interaction µ C in Eq. (46).(g) The averaged SF term µ s in Eq. (47). Note thatsince there is no computed data for Be and In with SOI,we are presenting data for these two elements withoutSOI. T c s of Zn and Cd (V, Nb, Tc, Pb) are overesti-mated (underestimated) in the differential plot (a) whilethose of Zn Cd, W, Ir (Be, V, Rh) are overestimated(underestimated) in the ratio plot (b). Therefore, forV, Zn, and Cd, the focus should be on examining theaccuracy of T c ; therefore, we attempted to find featuresof these three materials from elemental systems. FromFig. 8 (g), we can see that V has extremely large µ s .When the system has large SF, the SF-mediated inter-action in Eq. (33) changes rapidly because the inverse ofthe matrix in Eq. (48) approaches to a singular matrix.Therefore, the SF of such systems needs to computedmore precisely, for example, by including the gradientcollection into the magnetic exchange-correlation kernelin Eq. (25). However, it is difficult to identify the sig-nificant difference between Zn, Cd, and other materials.For example, the parameters of Ga are extremely closeto those of Zn. However, the calculated T c of Ga is ingood agreement with the experimental one.By leaving from the superconductivity, we can see thefollowing features from Fig. 8: DOS and µ C are show-ing very similar behavior. Since µ C in Eq. (46) can beapproximated to the DOS times the averaged Coulombinteraction, the synchronicity indicates that the elemen-tal materials in Fig. 8 have nearly the same screenedCoulomb repulsion. ω ln has peaks around group 6-9 onboth periods 5 and 6. Additionally, the frequencies atboth the peaks are extremely close, although the atomicmasses of these periods are different; these behaviors canbe explained by the following Friedel’s theory . The ma- terials in groups 6-9 have high cohesive energy because ofthe half-filled d-orbitals; the high cohesive energy leadsto the hardness and the high phonon frequencies of thematerials. Since 5d orbitals are more delocalized than4d orbitals, the binding energy increases for 5d materi-als. Therefore, the phonon frequency is unchanged be-cause of the cancellation between the atomic mass andthe stronger bonding. VI. SUMMARY
We performed the benchmark calculations of SCDFTusing our open-source software package
Superconduct-ing Toolkit , which uses our newly developed methodfor treating SOI together with SF. We presented bench-mark results of superconducting properties calculated bySCDFT for 35 elemental materials together with compu-tational details, and discussed the accuracy of the pre-dicted T c and the effects of SF and SOI up on T c . Wefound that the calculations, including SOI and SF, canquantitatively reproduce the experimental T c s. The SFis essential especially for the transition metals; still theeffect of SOI is small for elemental systems excepting Tc,Sn, Re, Tl, and Pb. We also reproduced the absenceof the superconductivity in the alkaline, alkaline earth,and noble metals. Our result could be used to check thevalidity of future work such as the high-throughput cal-culation for exploring new superconductors. Moreover,the knowledge of this benchmark can be used to improvethe methodology of SCDFT. For example, we can focuson Zn and Cd as a target for the next theoretical im-provement. It is straightforward to extend the currentbenchmark calculation into the binary, ternary, and qua-ternary systems. From such a benchmark, we check sys-tematically the accuracy of SCDFT for the compound su-perconductors such as magnesium diboride , cuprates ,iron-based , and heavy-fermion superconductors , etc.This type of benchmark calculation should be performedwhenever there is room for theoretical improvements sothat the applicability of SCDFT can be extended as auniversal tool. ACKNOWLEDGMENTS
We thank Ryosuke Akashi for the fruitful discussionabout the spin-fluctuation. This work was supported byPriority Issue (creation of new functional devices andhigh-performance materials to support next-generationindustries) to be tackled using Post ‘K’ Computer fromthe MEXT of Japan. Y. H. is supported by Japan So-ciety for the Promotion of Science through Program forLeading Graduate Schools (MERIT). The numerical cal-culations in this paper were done on the supercomputersin ISSP and Information Technology Center at the Uni-versity of Tokyo.4 ∗ [email protected] L. N. Oliveira, E. K. U. Gross, and W. Kohn, Phys. Rev.Lett. , 2430 (1988). M. L¨uders, M. A. L. Marques, N. N. Lathiotakis, A. Floris,G. Profeta, L. Fast, A. Continenza, S. Massidda, andE. K. U. Gross, Phys. Rev. B , 024545 (2005). E. R. Margine and F. Giustino, Phys. Rev. B , 024505(2013). W. L. McMillan, Phys. Rev. , 331 (1968). R. Dynes, Solid State Commun. , 615 (1972). A. Sanna, J. A. Flores-Livas, A. Davydov, G. Profeta,K. Dewhurst, S. Sharma, and E. K. U. Gross, Jour-nal of the Physical Society of Japan , 041012 (2018),https://doi.org/10.7566/JPSJ.87.041012. F. Essenberger, A. Sanna, A. Linscheid, F. Tandetzky,G. Profeta, P. Cudazzo, and E. K. U. Gross, Phys. Rev.B , 214504 (2014). M. A. L. Marques, M. L¨uders, N. N. Lathiotakis, G. Pro-feta, A. Floris, L. Fast, A. Continenza, E. K. U. Gross,and S. Massidda, Phys. Rev. B , 024546 (2005). A. Floris, G. Profeta, N. N. Lathiotakis, M. L¨uders,M. A. L. Marques, C. Franchini, E. K. U. Gross, A. Con-tinenza, and S. Massidda, Phys. Rev. Lett. , 037004(2005). A. Sanna, G. Profeta, A. Floris, A. Marini, E. K. U. Gross,and S. Massidda, Phys. Rev. B , 020511 (2007). G. Profeta, C. Franchini, N. N. Lathiotakis, A. Floris,A. Sanna, M. A. L. Marques, M. L¨uders, S. Massidda,E. K. U. Gross, and A. Continenza, Phys. Rev. Lett. ,047003 (2006). P. Cudazzo, G. Profeta, A. Sanna, A. Floris, A. Conti-nenza, S. Massidda, and E. K. U. Gross, Phys. Rev. Lett. , 257001 (2008). R. Akashi, M. Kawamura, S. Tsuneyuki, Y. Nomura, andR. Arita, Phys. Rev. B , 224513 (2015). F. Essenberger, A. Sanna, P. Buczek, A. Ernst, L. San-dratskii, and E. K. U. Gross, Phys. Rev. B , 014503(2016). R. Akashi and R. Arita, Phys. Rev. Lett. , 057006(2013). T. Nomoto, M. Kawamura, T. Koretsune, R. Arita,T. Machida, T. Hanaguri, M. Kriener, Y. Taguchi, andY. Tokura, Phys. Rev. B , 014505 (2020). A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, andI. Tanaka, Phys. Rev. Lett. , 205901 (2015). A. G¨orling and M. Levy, Phys. Rev. A , 196 (1994). R. Heid, K.-P. Bohnen, I. Y. Sklyadneva, and E. V.Chulkov, Phys. Rev. B , 174527 (2010). S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi,Rev. Mod. Phys. , 515 (2001). M. Gell-Mann and K. Brueckner, Phys. Rev. , 364(1957). E. K. U. Gross and W. Kohn, Phys. Rev. Lett. , 2850(1985). M. Kawamura, R. Akashi, and S. Tsuneyuki, Phys. Rev.B , 054506 (2017). E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984). J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. , 1175 (1957). J. Schrieffer,
Theory of Superconductivity , Advanced BookProgram Series (Advanced Book Program, Perseus Books,(1983)). P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B.Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli,M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso,S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti,A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerst-mann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. K¨u¸c¨ukbenli, M. Lazzeri, M. Marsili,N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O.de-la Roza, L. Paulatto, S. Ponc´e, D. Rocca, R. Saba-tini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov,I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu,and S. Baroni, Journal of Physics: Condensed Matter ,465901 (2017). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 1396 (1997). D. R. Hamann, Phys. Rev. B , 085117 (2013). M. Schlipf and F. Gygi, Computer Physics Communica-tions , 36 (2015). G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, andN. Marzari, npj Computational Materials , 72 (2018). M. Kawamura, Y. Gohda, and S. Tsuneyuki, Phys. Rev.B , 094515 (2014). http://sctk.osdn.jp/. J. Hamlin, Physica C: Superconductivity and its Applica-tions , 59 (2015), superconducting Materials: Conven-tional, Unconventional and Undetermined. K. A. Gschneidner, in
Solid State Physics: Advances inResearch and Applications , Vol. 16, edited by F. Seitz andD. Turnbull (Academic Press, 1964) p. 275. W. L. McMillan and J. M. Rowell, in
Superconductivity ,Vol. 1, edited by R. D. Parks (CRC Press, 1969) p. 561. K.-H. Lee, K. J. Chang, and M. L. Cohen, Phys. Rev. B , 1425 (1995). K.-H. Lee and K. J. Chang, Phys. Rev. B , 1419 (1996). A. Y. Liu and A. A. Quong, Phys. Rev. B , R7575(1996). S. Y. Savrasov and D. Y. Savrasov, Phys. Rev. B , 16487(1996). N. F. Berk and J. R. Schrieffer, Phys. Rev. Lett. , 433(1966). K. Tsutsumi, Y. Hizume, M. Kawamura, R. Akashi, andS. Tsuneyuki, To be submitted. J. Friedel, in
Theory of Magnetism in Transition Metals ,edited by Marshall, W. (Academic Press, 1967) p. 283. J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani,and J. Akimitsu, Nature (London) , 63 (2001). N. Plakida, in
High-Temperature Cuprate Superconductors:Experiment, Theory, and Applications , Springer Series inSolid-State Sciences, Vol. 166 (2010) pp. 1–570. G. R. Stewart, Rev. Mod. Phys. , 1589 (2011). C. Pfleiderer, Rev. Mod. Phys.81