Rigid muffin-tin approximation in plane-wave codes for fast modeling of phonon-mediated superconductors

TL;DR

This paper introduces a pseudopotential-based plane-wave implementation of the rigid muffin-tin approximation, enabling efficient evaluation of electron-phonon interactions for superconductors, suitable for high-throughput materials screening.

Contribution

It develops a computationally efficient pseudopotential-based RMTA method compatible with plane-wave codes, validated against full-potential calculations, facilitating rapid superconductor screening.

Findings

Excellent agreement with full-potential calculations

Enables extraction of atom- and symmetry-resolved electron-phonon parameters

Reduces computational cost for high-throughput superconductor discovery

Abstract

We present a pseudopotential-based plane-wave implementation of the rigid muffin-tin approximation (RMTA), offering a computationally efficient alternative to its traditional use in all-electron codes. This approach enables the evaluation of angular-momentum-resolved electron-phonon matrix elements and McMillan-Hopfield parameters of not only elemental transition metals but also their compounds. The results are benchmarked against full-potential linearized augmented plane wave calculations, showing excellent agreement. We further outline a practical route to extract atom- and symmetry-type-resolved electron-phonon coupling constants. By enabling the use of RMTA descriptors within high-throughput workflows, this framework significantly lowers the computational cost of screening candidate superconductors, providing a valuable tool for materials discovery.

Tables

Table 1. Table 1: Estimated semi-phenomenological EPC constant, λ \lambda , based on calculated McMillan-Hopfield parameters, η \eta , at specified MT-radii, r MT r_{\text{MT}} , and experimental Debye temperature, θ D \theta_{\text{D}} , for simple metals. θ D \theta_{\text{D}} values are taken from corresponding references. MH parameters η ref \eta_{\text{ref}} , calculated with the APW method and MJW potentials [ 14 ] , are provided for the reference.

Structure	$r_{\text{MT}}$	$\eta$	${\eta }_{\text{ref}}$ [14]	${\theta }_{\text{D}}$	$\lambda$
	(bohr)	$\left(\frac{\text{eV}}{{\text{Å}}^2}\right)$	$\left(\frac{\text{eV}}{{\text{Å}}^2}\right)$	(K)
V (BCC)	2.48	8.01	6.89	385 [73]	1.19
Cu (FCC)	2.41	0.57	0.46	343 [74]	0.09
Nb (BCC)	2.70	9.07	7.63	275 [75]	1.45
Mo (BCC)	2.57	5.56	5.80	430 [76]	0.35
Rh (FCC)	2.56	4.95	4.72	480 [77]	0.24
Pd (FCC)	2.66	2.65	1.99	299 [78]	0.31
Ag (FCC)	2.73	0.28	0.29	225 [74]	0.06

Equations

$$\Delta V_{\kappa }(r)\approx -\Delta {\mathbf{\tau }}_{\kappa }\cdot \nabla V_{\kappa }(r)=-\left(\Delta {\mathbf{\tau }}_{\kappa }\cdot \hat{\mathbf{r}}\right)\frac{d{V}_{\kappa }(r)}{dr}.$$

$$\lambda$$

$${\tilde{\lambda }}_{\kappa }=\frac{{\tilde{\eta }}_{\kappa }}{M_{\kappa }}⟨\frac{1}{{\omega }_{\kappa }^2}⟩\approx \frac{{\tilde{\eta }}_{\kappa }}{M_{\kappa }⟨{\omega }_{\kappa }^2⟩}.$$

$${\tilde{\eta }}_{\kappa }$$

$$\times {\int }_{\text{BZ}}d{\mathbf{k}}^{\prime }\delta \left({\epsilon }_{j,{\mathbf{k}}^{\prime }}-{\epsilon }_{\text{F}}\right)$$

$$\times {\left|⟨{\Psi }_{j,{\mathbf{k}}^{\prime }}|\nabla V_{\kappa }(r)|{\Psi }_{i,\mathbf{k}}⟩\right|}^2,$$

$${\lambda }_{\mu }$$

$$\lambda$$

$${\eta }_{\mu }=N_{\text{atoms}}^{\mu }{\tilde{\eta }}_{\kappa }\forall \kappa \in \{\mu \}.$$

$${\Psi }_{i,\mathbf{k}}(\mathbf{r})$$

$$a_{i\kappa ,lm}(\mathbf{k})=a_{i\kappa ,l}(k)Y_{lm}^{\ast }(\hat{\mathbf{k}}),$$

$${\tilde{\eta }}_{\kappa }$$

$$\times \frac{n_{\kappa ,l}(r_{\kappa ,\text{MT}},{\epsilon }_{\text{F}})\times n_{\kappa ,l+1}(r_{\kappa ,\text{MT}},{\epsilon }_{\text{F}})}{N({\epsilon }_{\text{F}})},$$

$$M_{\kappa ;l,l+1}(r_{\kappa ,\text{MT}},{\epsilon }_{\text{F}})=$$

$$\frac{{\int }_0^{r_{\kappa ,\text{MT}}}dr{r}^2{R}_{\kappa ,l}(r,{\epsilon }_{\text{F}})\frac{d{V}_{\kappa }(r)}{dr}{R}_{\kappa ,l+1}(r,{\epsilon }_{\text{F}})}{\sqrt{{\int }_0^{r_{\kappa ,\text{MT}}}dr{r}^2{R}_{\kappa ,l}^2(r,{\epsilon }_{\text{F}}){\int }_0^{r_{\kappa ,\text{MT}}}dr{r}^2{R}_{\kappa ,l+1}^2(r,{\epsilon }_{\text{F}})}}.$$

$$n_{\kappa ,l}$$

$$=\sum _m\frac{\Omega }{(2\pi {)}^3}{\int }_0^{r_{\kappa ,\text{MT}}}dr{r}^2{R}_{\kappa ,l}^2(r,{\epsilon }_{\text{F}})$$

$$\times \sum _i{\int }_{\text{BZ}}d\mathbf{k}\delta ({\epsilon }_{i,\mathbf{k}}-{\epsilon }_{\text{F}}){\left|a_{i\kappa ,lm}(\mathbf{k})\right|}^2,$$

$$M_{\kappa ;l,l+1}(r_{\text{MT}},{\epsilon }_{\text{F}})=$$

$$\displaystyle\frac{[V_{\kappa}(r)-\varepsilon_{\text{F}}]r^{2}-(L_{\kappa,l}-l)[L_{\kappa,l+1}+(l+2)]}{r\sqrt{\frac{\partial L_{\kappa,l}}{\partial\varepsilon}\cdot\frac{\partial L_{\kappa,l+1}}{\partial\varepsilon}}}\Bigg{|}_{\begin{subarray}{c}\varepsilon=\varepsilon_{\text{F}}\\ r=r_{\kappa,\text{MT}}\end{subarray}},$$

$$L_{\kappa ,l}(r,\epsilon )=r\frac{\left(\frac{\partial R_{\kappa ,l}(r,\epsilon )}{\partial r}\right)}{R_{\kappa ,l}(r,\epsilon )}.$$

$$\tilde{V}(\mathbf{G})$$

$$V_{\kappa ,\text{loc}}(r)$$

$$=\sum _{\mathbf{G}}{e}^{i\mathbf{G}\cdot {\mathbf{\tau }}_{\kappa }}\tilde{V}(\mathbf{G})\frac{\sin Gr}{Gr}.$$

$$\displaystyle\left[-\frac{d^{2}}{dr^{2}}+V^{l}_{\kappa,\text{eff}}(r)\right]\left[rR_{\kappa,l}(r,\varepsilon)\right]=\varepsilon\left[rR_{\kappa,l}(r,\varepsilon)\right]\Big{|}_{\varepsilon=\varepsilon_{\text{F}}},$$

$$V_{\kappa ,\text{eff}}^l(r)$$

$$=\left[V_{\kappa ,\text{loc}}(r)+V_{\kappa ,\text{SL}}^l(r)\right]+\frac{l(l+1)}{r^2}.$$

$$\begin{cases}\Psi_{i,\bm{k}}(\bm{r})&=\sum_{lm}a_{\kappa i,lm}(\bm{k})R_{\kappa,l}(r_{\kappa,\text{MT}},\varepsilon_{i,\bm{k}})Y_{lm}(\hat{\bm{r}}),\\ \Psi_{i,\bm{k}}(\bm{r})&=\psi_{i,\bm{k}}(\bm{\tau}_{\kappa}+\bm{r})\Big{|}_{r=r_{\kappa,\text{MT}}}\\ &=\frac{1}{\sqrt{\Omega}}\sum_{\bm{G}}e^{i(\bm{k}+\bm{G})\cdot\left(\bm{\tau}_{\kappa}+\bm{r}\right)}\tilde{\psi}_{i}(\bm{k}+\bm{G})\Bigg{|}_{r=r_{\kappa,\text{MT}}},\end{cases}$$

$$a$$

$$\displaystyle\times\int d\hat{\bm{r}}\left[\sum_{\bm{G}}e^{i(\bm{k}+\bm{G})\cdot\left(\bm{\tau}_{\kappa}+\bm{r}\right)}\tilde{\psi}_{i}(\bm{k}+\bm{G})\right]Y^{\ast}_{lm}({\hat{\bm{r}}})\Bigg{|}_{r=r_{\kappa,\text{MT}}}.$$

$$n_{\kappa ,l}(r_{\kappa ,\text{MT}},{\epsilon }_{\text{F}})=\sum _{m=-l}^l{n}_{\kappa ,lm}(r_{\kappa ,\text{MT}},{\epsilon }_{\text{F}})$$

$$=\sum _{m=-l}^l\frac{r_{\kappa ,\text{MT}}}{(2\pi {)}^3}\left(-\frac{\partial }{\partial \epsilon }{L}_{\kappa ,l}\right){\int }_{\text{BZ}}d\mathbf{k}\sum _{i=1}^{N_{\text{states}}}\delta ({\epsilon }_{i,\mathbf{k}}-{\epsilon }_{\text{F}})\times$$

$$\displaystyle\left|\int d\hat{\bm{r}}\left[\sum_{\bm{G}}e^{i(\bm{k}+\bm{G})\cdot(\bm{\tau}_{\kappa}+\bm{r})}\tilde{\psi}_{i}(\bm{k}+\bm{G})\right]Y^{\ast}_{lm}({\hat{\bm{r}}})\right|^{2}\Bigg{|}_{\begin{subarray}{c}r=r_{\kappa,\text{MT}}\\ \varepsilon=\varepsilon_{\text{F}}\end{subarray}}.$$

$$N({\epsilon }_{\text{F}})=\frac{\Omega }{(2\pi {)}^3}{\int }_{\text{BZ}}d\mathbf{k}\sum _i\delta \left({\epsilon }_{i,\mathbf{k}}-{\epsilon }_{\text{F}}\right).$$

$${\tilde{\eta }}_{\kappa }$$

$$\times {\int }_{\text{BZ}}d{\mathbf{k}}^{\prime }\delta \left({\epsilon }_{j,{\mathbf{k}}^{\prime }}-{\epsilon }_{\text{F}}\right){\int }_{\text{MT}}d\mathbf{r}{\int }_{\text{MT}}d{\mathbf{r}}^{\prime }$$

$$\times {\Psi }_{j,{\mathbf{k}}^{\prime }}({\mathbf{r}}^{\prime }){\Psi }_{i,\mathbf{k}}^{\ast }({\mathbf{r}}^{\prime }){\Psi }_{j,{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r}){\Psi }_{i,\mathbf{k}}(\mathbf{r})$$

$$\times \left(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime }\right)\frac{d{V}_{\kappa }(r)}{dr}\frac{d{V}_{\kappa }(r^{\prime })}{d{r}^{\prime }}.$$

Full text

I Introduction

Historically, the first practical method for calculating the electron-phonon coupling (EPC) constant in metals was rigid muffin-tin approximation (RMTA) [1, 2, 3], in which the electronic potential within a sphere around each atom, called the muffin-tin (MT) sphere, was assumed to follow ionic displacements rigidly. Since the RMTA was formulated in terms of MT potentials and wavefunctions expanded in radial functions and spherical harmonics, it was naturally implemented in many all-electron methods, such as the augmented plane wave (APW) [4, 5, 6, 7, 8], Korringa-Kohn-Rostoker (KKR) [9, 10, 11], and linear muffin-tin orbital (LMTO) [12, 13, 6] methods, where atomic functions inside the MT sphere are evaluated explicitly. The RMTA method worked reasonably well for materials with localized $d$-orbitals and close-packed structures (e.g., V, Nb, Mo, and Pd) [12, 14, 15, 16, 17, 18, 19, 20], and was actively used to study superconductivity of various materials, such as $A15$ [21, 16] and $C15$ [22] compounds, hydrides [23], nitrides [24], etc. While this method was designed with transition metals in mind, and routinely tended to underestimate the EPC in $sp$-metals, corrections accounting for the long-range tails of the ionic MT potential at large distances in such materials have been proposed [25, 19] and successfully applied to, e.g., Ag [19] and Al, Ga, and Pb [25].

Later on, the RMTA was superseded by the much more accurate and universal linear-response method, formulated within density-functional perturbation theory (DFPT), [26, 27, 28, 29] and by superconducting density-functional theory (SCDFT) [30, 31, 32, 33], both of which are well-suited for fast pseudopotential-based plane-wave (PS$+$PW) codes, and are widely used for predicting superconducting properties. The first high-throughput materials screenings involving hundreds of $\alpha^{2}F(\omega)$-based superconductivity calculations have been recently reported [34, 35, 36, 37], but the high DFPT computational cost imposes limitations in the form of restricted unit cell sizes or confined composition spaces. While these efforts mark important milestones in applying machine learning to guide superconductor discovery, larger datasets, producible with fast physics-based methods, are needed for identifying broader design rules.

Several fast descriptors for superconductors have been proposed, typically tailored to particular materials families. For MgB2-type layered compounds, the difference between in-plane boron phonon frequencies at the M and $\Gamma$ points was used as a softening marker expected to correlate with the EPC [38]. In $sp$-metals, scaled EPC constants evaluated at the Brillouin-zone center were introduced as a proxy for the total EPC [39, 40], but this approach cannot describe cases where zone-boundary phonons play a significant role. More recently, frozen-phonon calculations at a small number of high-symmetry $q$-points combined with energy band shifts were proposed to estimate the EPC across materials of general chemistry, with particular emphasis on phases with large unit cells [41], although to date the method has been applied exclusively to $sp$-bonded metals and hydrides. Real-space descriptors have also been explored, including the electron-localization function [42, 43] and spatial derivatives of the Kohn–Sham potential [44], which were used to estimate $T_{c}$ in hydrogen-based superconductors. While effective within their target classes, usually comprised of light elements, these descriptors underscore the need for approaches that can extend to transition metals and their compounds.

The RMTA could address this problem since it has been successfully applied to transition metals, but so far it has not been formulated in a way suitable for pseudopotential-based codes, which are preferably used by modern search algorithms. To assist future superconductor predictions and enhance the functionality of such codes, the current work presents a practical approach for evaluating EPC constants using the RMTA within the plane-wave pseudopotential formalism. Developed method, which is a fast postprocessing step on top of a standard self-consistent field (SCF) density functional theory (DFT) calculation, is applied to transition metals and their compounds. Calculated electron-phonon (EP) matrix elements [12] and McMillan–Hopfield (MH) parameters [45, 46, 3, 12, 17] are shown to agree well with the RMTA implementation within the full-potential linearized augmented plane wave (FLAPW) method.

The paper is organized as follows: Sec. II introduces the notations and describes the methodology, with more detailed derivations provided in the Appendices. Sec. III outlines relevant technical aspects of the tests. Calculated partial EP matrix elements and MH parameters are compared to those obtained with an FLAPW code in Sec. IV. In the same section, estimates of EPC constants for simple metals are obtained using experimental Debye temperatures. Sec. V discusses how the RMTA descriptors can be leveraged with machine learning algorithms. Unless stated otherwise, atomic Rydberg units are assumed throughout the paper.

II Methodology

This section describes three underlying assumptions involved in the RMTA [47, 48] along with subtleties of atom- and type-resolved EPC constants and MH parameters, generalizes standard RMTA equations to the case of multiple symmetry types, introduces partial EP matrix elements and partial electronic density of states (eDOS) needed to calculate MH parameters, and summarizes key equations derived in this work that make the RMTA compatible with pseudopotentials.

II.1 Assumption 1: Rigid-ion approximation

The RMTA relies on the assumption that a change in the potential acting on an electron located at the vector $\bm{r}$ from the ionic center $\bm{\tau}_{\kappa}$ is given by the rigid displacement of the electron-ion interaction potential $V_{\kappa}(\bm{r})$ (see Fig. 1), given by the dot product of the ionic displacement $\Delta\bm{\tau}_{\kappa}$ and gradient of the potential $\bm{\nabla}V_{\kappa}$. Using the spherical part of the potential, $V_{\kappa}(r)$, present work evaluates this change in Eq. (1), as was done in the original RMTA formulations [1, 2, 12].

[TABLE]

II.2 Assumption 2: Local-vibration approximation

In the RMTA, the EP interaction is treated as local in real space [1, 3], such that the total EPC constant for the whole material, $\lambda$, can be represented as a sum of atomic EPC constants $\tilde{\lambda}_{\kappa}$ over $N_{\text{atoms}}$ atoms in the considered cell [49, 50]:

[TABLE]

The EPC constant for each atom is expressed as a ratio [see Eq. (3)] of an “electronic” contribution, represented by MH parameter $\tilde{\eta}_{\kappa}$, and a “phonon” part, given by the mass of the ion $M_{\kappa}$ times its average phonon frequency $\left<\omega_{\kappa}^{2}\right>$ [45, 46, 3, 17].

[TABLE]

Although there are known approaches [3, 51, 17] to avoid full phonon calculations to obtain $\left<\omega_{\kappa}^{2}\right>$, their implementation is out of the scope of the present work. Instead, as an example, we approximate $\left<\omega_{\kappa}^{2}\right>$ in simple metals using experimental Debye temperatures $\theta_{\text{D}}$, via the relation $\left<\omega_{\kappa}^{2}\right>\approx\frac{1}{2}\theta_{\text{D}}^{2}$.

The main focus of this work is the derivation of the expressions for the “electronic” part, which are more analytically involved than those for the “phonon” part, as they directly depend on pseudopotentials and atomic functions derived from them. The corresponding atomic MH parameter, $\tilde{\eta}_{\kappa}$, is defined as the Fermi-surface average of the squared EP matrix elements for states at the Fermi level, divided by the eDOS at the Fermi energy:

[TABLE]

where $N(\varepsilon_{\text{F}})$ is the eDOS at the Fermi level per cell per spin; $\varepsilon_{\text{F}}$ is Fermi energy; $\varepsilon_{i,\bm{k}}$ and $\varepsilon_{j,\bm{k}^{\prime}}$ are the Hamiltonian eigenvalues of the states $i$ and $j$ corresponding to the wavefunctions $\Psi_{i,\bm{k}}(\bm{r})$ and $\Psi_{j,\bm{k}^{\prime}}(\bm{r})$ at momenta $\bm{k}$ and $\bm{k}^{\prime}$, respectively; $N_{\text{states}}$ is the total number of electronic states; $\Omega$ is the unit cell volume; and “BZ” marks the integration over the Brillouin zone. The eDOS in the denominator of Eq. (4) makes the atomic MH parameter and the corresponding atomic EPC constant dependent on the cell size. The total EPC constant, $\lambda$ in Eq. (2), however, is independent of the cell size.

In materials where different Wyckoff positions are occupied, it is convenient to introduce the EPC constant $\lambda_{\mu}$ for each of the $N_{\text{types}}$ symmetry types, as defined in Eq. (LABEL:eq:lambda_mu). Each symmetry type includes only atoms of the same chemical element that are related by crystallographic symmetry operations, with $N_{\text{atoms}}^{\mu}$ being the number of atoms per type $\mu$. The second equality of Eq. (LABEL:eq:lambda_mu) follows from the fact that all atoms belonging to the same type $\mu$ share identical values of $\tilde{\lambda}_{\kappa}$.

[TABLE]

Unlike the atomic EPC constants, the type-specific EPC constants, $\lambda_{\mu}$, are independent of the cell size. The total EPC constant is then given by a sum over the EPC constants corresponding to different symmetry types [see Eq. (6)] [16, 52, 53, 54, 55]. Importantly, the ability to find the contributions of individual symmetry types to the overall EPC in compounds allows for a straightforward comparison between these contributions in different materials.

[TABLE]

Since the masses $M_{\kappa}$ and the averaged phonon frequencies $\left<\omega_{\kappa}^{2}\right>$ are identical for atoms of the same symmetry type, MH parameters for different symmetry types can also be introduced as

[TABLE]

Similar to the type-specific EPC constants, the type-specific MH parameters do not depend on the cell size. Although the present formalism uses atomic parameters, only the type-specific ones are reported in Section IV since they can be compared across different materials regardless of the unit cell choice.

II.3 Assumption 3: Spherical-band approximation

The RMTA assumes that the wavefunction for the state $i$ inside the MT sphere of atom $\kappa$, with radius $r_{\kappa,\text{MT}}$, can be expanded in terms of radial functions $R_{\kappa,l}(r,\varepsilon_{i,\bm{k}})$ corresponding to the local screened potential $V_{\kappa}(r)$ and spherical harmonics $Y_{lm}(\hat{\bm{r}})$ characterized by quantum numbers $l$ and $m$, with momentum-dependent coefficients $a_{i\kappa,lm}(\bm{k})$:

[TABLE]

As shown in Appendix A, the spherical-band approximation that separates the radial ($k$) and angular ($\hat{\bm{k}}$) dependencies on the wavevector $\bm{k}$ as

[TABLE]

along with the separation of radial ($r$) and angular ($\hat{\bm{r}}$) components of $\bm{r}$, is crucial. This approximation enables a significant simplification of the MH parameter $\tilde{\eta}_{\kappa}$, reducing Eq. (4) to

[TABLE]

where the partial EP matrix elements are defined as an integral of two normalized radial functions with the potential gradient over the MT-sphere volume,

[TABLE]

As a result of the spherical-band approximation, Eq. (10) includes only partial EP matrix elements corresponding to dipole transitions with $\Delta l=\pm 1$. The quantities $n_{\kappa,l}$ represent the partial eDOS for angular momentum $l$ inside the MT sphere of atom $\kappa$, defined as

[TABLE]

It is important to point out that the spherical-band approximation is not required for the partial eDOS calculation, as can be seen in Eq. (12), where general coefficients $a_{i\kappa,lm}(\bm{k})$ are used.

Hereafter, the task of computing MH parameters in PS+PW codes comes down to the calculation of partial EP matrix elements and partial eDOS.

II.4 Partial EP matrix elements evaluated

on the MT-sphere

The original expression for the partial EP matrix elements, introduced by Gaspari and Gyorffy [Eq. (11)] [1, 2], involves integration of the radial functions $R_{\kappa,l}(r,\varepsilon_{\text{F}})$ and $R_{\kappa,l+1}(r,\varepsilon_{\text{F}})$, along with the derivative of the MT potential $dV_{\kappa}(r)/dr$, over the volume of the MT-sphere [Fig. 2(a)]. This expression has three problems for the PS+PW codes: (1) A pseudopotential (PS) differs from the all-electron (AE) potential within a cut-off radius $r_{\text{c}}$, where pseudization takes place, rendering the integration of the PS derivative in that region unreliable; (2) The screened self-consistent potential, as calculated in PS+PW codes, is periodic in the real space, and needs to be converted to localized screened atomic potentials; (3) The radial functions $R_{\kappa,l}(r,\varepsilon)$ are not calculated explicitly and need to be derived from available DFT data.

In our method, the first problem is solved by referring to the following reformulated expression for Eq. (11):

[TABLE]

It was introduced by Pettifor [12] based on the radial Schrödinger equation for radial functions $R_{\kappa,l}(r,\varepsilon)$ and their normalization conditions, provided in Appendix B. Here, $V_{\kappa}(r)$ is the local screened potential of the atom $\kappa$, $L_{\kappa,l}(r,\varepsilon)$ is the logarithmic derivative of the radial part of the wavefunction evaluated around the Fermi level, and $\frac{\partial L_{\kappa,l+1}}{\partial\varepsilon}$ is the energy derivative of the logarithmic derivative. In this expression, the values of the potential and logarithmic derivatives at the MT-radius are sufficient to calculate the partial EP matrix elements.

The second problem is addressed by using properties of the Fourier transformation. The self-consistent DFT pseudopotential is given as a periodic local screened potential $V_{\text{loc}}(\bm{r}_{0})$ [Fig. 3(b)] on the 3D real-space grid $\bm{r}_{0}$, with a center at an arbitrary origin of the calculated unit cell. The Fourier transform of this potential is represented by coefficients

[TABLE]

where $\bm{G}$ is a reciprocal lattice vector.

A spherically-symmetric local screened atomic potential, $V_{\kappa,\text{loc}}(r)$, is then expressed in Eq. (15) by obtaining coefficients of the plane-wave expansion of the periodic potential on a shifted grid, $\bm{r}$, that is centered on a specific atom at $\bm{\tau}_{\kappa}$ [Fig. 3(a-b)], as well as by using the expansion of $e^{i\bm{G}\cdot\bm{r}}$ in spherical Bessel functions and spherical harmonics for only $l=m=0$ [see Appendix C for details].

[TABLE]

In this equation, $r~=~|\bm{r}|$, $G~=~|\bm{G}|$, and $e^{i\bm{G}\cdot\bm{\tau}_{\kappa}}\tilde{V}(\bm{G})$ are modified coefficients of the plane-wave expansion of the same periodic potential $V_{\text{loc}}(\bm{r})$ [Fig. 3(b)] but on a shifted radial grid, $\bm{r}$, that is centered on the atom. The potential $V_{\text{loc}}(\bm{r})$ is DFT-screened. Hence, the extracted local pseudopotential $V_{\kappa,\text{loc}}(r)$ is screened, too. Note that the local PS is equal to the AE potential above the cut-off radius: $V_{\kappa}(r)=V_{\kappa,\text{loc}}(r)~\forall~r>r_{\text{c}}$.

For the third problem, the method of choice is to generate radial functions $R_{\kappa,l}(r,\varepsilon)$, which are compatible with the extracted screened potential $V_{\kappa,\text{loc}}(r)$, for a particular energy $\varepsilon$ by solving the corresponding radial Schrödinger equation,

[TABLE]

with an effective potential $V^{l}_{\kappa,\text{eff}}(r)$ containing $V_{\kappa,\text{loc}}(r)$, $l$-dependent semilocal (SL) potential $V^{l}_{\kappa,\text{SL}}(r)$ that is non-zero only in the core region $r<r_{\text{c}}$ (see Fig. 4), and centrifugal potential $l(l+1)/r^{2}$:

[TABLE]

Eq. (16) is solved numerically at the specified energy $\varepsilon$ by outward integration from $r$ value close to zero to the desired MT-radius. This way, the energy derivatives of the logarithmic derivatives $\partial L_{\kappa,l}/\partial\varepsilon$ can be obtained by solving the radial Schrödinger equation at different energies close to the Fermi level. The contribution from the SL potential is pertinent to the PS approach: Different pseudopotentials are needed to generate radial functions for different $l$ values, while only one AE potential would generate all radial functions.

Although most PS+PW codes use pseudopotentials in the non-local (NL) (dependent on $r$ and $r^{\prime}$) Vanderbilt-Kleinman-Bylander (VKB) form [56, 57, 58] due to its computational efficiency in reciprocal space, we use the SL form for simplicity in the current implementation. For the non-local ONCVPSP pseudopotentials by Hamann [59, 60] utilized in this work, automatic conversion of the NL components, read directly from the pseudopotential file, to the SL form is straightforward and described in Appendix D. In principle, solving the radial Schrödinger equation with NL potentials directly [61] would give the same results.

II.5 Partial eDOS inside the MT-sphere and total eDOS

The partial eDOS inside the MT-sphere of atom $\kappa$ are derived based on the assumption that on the surface of the MT-sphere the wavefunction of the state $i$ can be expanded in both spherical harmonics and plane waves at the same time:

[TABLE]

where $\tilde{\psi}_{i}(\bm{k}+\bm{G})$ are coefficients of the Fourier expansion of the wavefunction $\Psi_{i,\bm{k}}(\bm{r})$.

From the system of Eqs. (18), the coefficients $a_{\kappa i,lm}(\bm{k})$ are found following Eq. (19), which is calculated by using the Gauss-Legendre quadrature method [62].

[TABLE]

Then, using Eqs. (19) and (25) (see Appendix A), the partial eDOS, $n_{\kappa,l}$, can be evaluated from the Fourier coefficients, $\tilde{\psi}_{i}(\bm{k}+\bm{G})$, and the logarithmic derivatives, $\frac{\partial}{\partial\varepsilon}L_{\kappa,l}$, by Eq. (20). Note that this method allows for the calculation of not only the $l$-resolved eDOS, but also $m$-resolved eDOS, $n_{\kappa,lm}$.

[TABLE]

As can be seen from Eq. (13), the squared partial EP matrix element, $M^{2}_{l,l+1}$, contains energy derivatives of the logarithmic derivatives, $\frac{\partial}{\partial\varepsilon}L_{\kappa,l}$ and $\frac{\partial}{\partial\varepsilon}L_{\kappa,l+1}$, in the denominator. These energy derivatives are normalization integrals for the radial functions $R_{l}$ and $R_{l+1}$ (see Appendix B) that cancel out when $M^{2}_{l,l+1}$ is multiplied by the partial eDOS $n_{l}$ and $n_{l+1}$ from Eq. (20). Hence, if the calculation of the partial eDOS themselves is not required, Eq. (10) for the MH parameters can be rewritten in a form independent of $\frac{\partial}{\partial\varepsilon}L_{\kappa,l}$ and $\frac{\partial}{\partial\varepsilon}L_{\kappa,l+1}$, such that energy derivatives of the logarithmic derivatives, or normalization integrals, do not have to be computed.

The eDOS at the Fermi level per cell per spin is computed as a Brillouin-zone integral of the delta-functions summed over all electronic states,

[TABLE]

III Technical details

The present implementation was developed within the Quantum ESPRESSO (QE) package [63, 64, 65], which was employed in this study for calculating the self-consistent pseudopotentials with the plane-wave (PW) code. Norm-conserving ONCVPSP [59, 60] pseudopotentials with the Perdew-Burke-Ernzerhof (PBE) [66] generalized-gradient approximation (GGA) for the exchange-correlation functional were used for all atoms, with a wavefunction energy cutoff of 120 Ry. The RMTA code was executed as a postprocessing step after the self-consistent field (SCF) calculation. For the quantitative verification presented in Sec. IV, all structures were modeled with their primitive unit cells, and MT-radii were chosen to correspond to touching spheres. The default three-dimensional FFT mesh for charge density and self-consistent potential, automatically determined by the energy cutoff, was retained. A $24\times 24\times 24$ Monkhorst-Pack $\bm{k}$-point grid and Methfessel-Paxton smearing [67] with $\sigma=10^{-2}$ Ry were utilized for the approximation of the occupation function in the SCF calculation. The tetrahedron method of integration [68] on the same $\bm{k}$-point grid was selected for the evaluation of the partial and total eDOS across all structures, since it showed faster convergence with respect to the $\bm{k}$-point grid than smeared delta-function approximations (see Section II of the Supplementary Information [69]).

Calculated partial EP matrix elements and McMillan-Hopfield parameters were compared to the FLAPW code flair [70]. In FLAPW, semicore states of all atoms were treated as core states and included explicitly in core-orthogonalization. Equations (10) and (13) were evaluated with the atomic potential, radial functions, and partial and total eDOS calculated using the FLAPW method. The same MT-radii, $\bm{k}$-point grid, and smearing parameters as those used for the RMTA calculations with the PW code were employed to ensure consistency.

IV Quantitative verification

A key difference in how the partial EP matrix elements and MH parameters depend on the choice of the MT-radius is illustrated in Figs. 5(a)-(b) for the BCC structure of Nb. For all $l$ channels, $M^{2}_{l,l+1}$ decrease as MT-radius increases [Fig. 5(a)], consistent with Pettifor’s observation [12] that partial EP matrix elements alone are not sufficient to describe the electronic part of the EPC constant because they are strongly dependent on the MT-radius.

On the contrary, the partial MH parameters, which are presented in eV / Å2 following previous studies [12, 14, 16], appear to be less sensitive [Fig. 5(b)] to the choice of the MT-radius since they include the partial eDOS, which increases as the MT sphere expands. Notably, while the $df$ channel may seem dominant based on its corresponding partial EP matrix element, it contributes roughly the same to the total MH parameter $\eta$ as the $pd$ channel at the reference MT radius. Achieving a nearly-flat behavior of the total atomic MH parameter around the reference MT-radius was one of the objectives of this work. Nevertheless, the slight dependence on $r_{\text{MT}}$ indicates that careful selection of the MT-radius based on the cell geometry is advisable.

The good agreement between the PS+PW and FLAPW methods for the electronic part calculated in selected simple metals and compounds is demonstrated in Figs. 6 and 7, where Strukturbericht symbols [71] are used as space-group identifiers for compounds. Crystallographic parameters of the structures are provided in the Supplementary Information [69]. The differences in partial EP matrix elements ($pd$ and $df$) for metal atoms in the selected structures do not exceed 5% [Fig. 6], with PS+PW and FLAPW results almost overlapping. In Fig. 7, which shows MH parameters for symmetry types of the metal atoms, both partial and total MH parameters exhibit similar trends across the two methods. Both methods demonstrate that $\eta$ parameter of Nb atoms in NbN is mostly dominated by the $df$ channel contribution, in agreement with previous studies [16, 21]. The observed discrepancies bewteen PS+PW and FLAPW are primarily due to slight differences in the electronic density of states. Importantly, there is no tendency for one of the methods to produce consistently higher values than the other. Comparison of the total MH parameters, shown as circles in Fig. 7, to some previous APW results [14, 16] is discussed in Section III of the Supplementary Information [69].

Table 1 compares the calculated MH parameters for simple metals with the results of Papaconstantopoulos et al. [14], and uses them to estimate the EPC constant by approximating the “phonon” contribution via the experimental Debye temperature. It should be noted that the MH parameters from Ref. [14] were obtained with the augmented plane wave (APW) method and Moruzzi-Janak-Williams (MJW) LDA potentials [72]. The noticeable differences from the present results are likely due to the evolution of potentials and the use of different exchange-correlation functionals. While the prior MH parameters tend to be lower than those obtained in this work, the overall trend across elemental metals remains similar. This consistency suggests that the RMTA can be effectively used to identify materials with a higher likelihood of superconductivity, provided that the same RMTA implementation is applied consistently across all candidates.

V Discussions

Machine learning studies [54, 79, 80] have demonstrated that McMillan-Hopfield parameters and partial eDOS can serve as effective descriptors for identifying potential superconductors. In this context, the present method is well-suited for high-throughput screening of candidate materials. As a minimal criterion, materials with very low MH parameters can be reliably excluded from consideration, as they are unlikely to exhibit superconductivity.

For a more complete estimate of the EPC constant, the phonon part, given by the averaged squared phonon frequencies $\left<\omega^{2}_{\kappa}\right>$, still needs to be evaluated using a computationally efficient approach. One practical strategy is to replace the direct evaluation of $\left<\omega^{2}_{\kappa}\right>$ with that of the interatomic force constants. These can be approximated from the diagonal elements of the inverted dynamical matrix, obtained at a few high-symmetry points on the Brillouin zone boundary via phonon linear-response calculations [3, 41], or from the diagonal elements of the inverse force-constant matrix at the $\Gamma$-point in frozen-phonon calculation, assuming a sufficiently large supercell [51, 17, 47].

To improve RMTA predictions for $sp$-metals, corrections to the asymptotic behavior of the local ionic potential could be introduced, as done by Zdetsis et al. [25] and Mazin et al. [19].

VI Conclusions

We presented a method for implementing the rigid muffin-tin approximation in pseudopotential-based plane-wave (PS+PW) codes to calculate the electronic part of the EPC constant, expressed through the McMillan-Hopfield parameter. The calculated values for a variety of simple metals and compounds were validated against results from a full-potential augmented plane wave code. The excellent agreement in partial EP matrix elements and McMillan-Hopfield parameters between the pseudopotential and all-electron methods confirms the reliability of the developed methodology. Some discrepancies with older references likely reflect the evolution of potentials and exchange-correlation functionals over time. As highlighted in previous machine learning studies, McMillan-Hopfield parameters computed within the PS+PW framework can be used as descriptors in machine learning algorithms for high-throughput screening of potential superconducting materials.

Acknowledgements.

This work was supported by National Science Foundation (Grants No. DMR-2320074, DMR-2320073, and DMREF-2323857).

Appendix A Assumptions of spherical potential and bands for McMillan-Hopfield parameters

This section demonstrates how the assumptions of a spherical potential and spherical bands lead to simplified expressions for the MH parameters. Following [1, 2, 12], Eq. (4) can be rewritten as

[TABLE]

While the explicit use of the wavefunction expansion from Eq. (8) would lead to four distinct sets of $\{l,m\}$ quantum numbers, corresponding to each wavefunction, the adoption of the spherical approximation in Eq. (9) and the orthonormality of the $\hat{\bm{k}}$- and $\hat{\bm{k}}^{\prime}$-dependent spherical harmonics reduce the number of independent angular momentum indices to two sets: $\{l,m\}$ and $\{l^{\prime},m^{\prime}\}$. As a result, Eq. (22) simplifies to

[TABLE]

The next simplification is achieved by taking advantage of the separation of the spatial dependence of $\bm{r}$ into radial, $r$, and angular, $\hat{\bm{r}}$, components. Applying the identity [81]

[TABLE]

along with the properties of Gaunt coefficients (integrals involving three spherical harmonics [82]), Eq. (23) simplifies to a form involving only dipole transitions with $l^{\prime}=l\pm 1$ and $m^{\prime}=m+m^{\prime\prime}$. The resulting expression for $\tilde{\eta}_{\kappa}$ is given by Eqs. (10), (11), and (12).

Notably, the partial eDOS, defined in Eq. (11), can be further simplified by using the normalization of the radial functions from Appendix B:

[TABLE]

In Eq. (25), the integrals over $r$ are replaced by energy derivatives of the logarithmic derivatives, ${\partial L_{\kappa,l}}/{\partial\varepsilon}$, evaluated only on the surface of the MT-sphere.

Appendix B Normalization of the radial functions

From the radial Schrödinger equation and its energy derivative, an integral of the squared radial function can be converted into the energy-derivative of the logarithmic derivative of this function as

[TABLE]

Appendix C Local screened atomic potential

The coefficients $\tilde{V}(\bm{G})$ in the plane-wave expansion of the periodic local screened potential $V_{\text{loc}}(\bm{r}_{0})$ [Fig. 3(b)] are given by Eq. (14), where the vector $\bm{r}_{0}$ starts at the origin of the calculated periodic unit cell [see Fig. 3(a)]. To obtain the periodic local potential on the atom-centered $\bm{r}$-grid, with the atom located at $\bm{\tau}_{\kappa}$, such that $\bm{r}_{0}=\bm{\tau}_{\kappa}+\bm{r}$, one can apply the inverse Fourier transform

[TABLE]

In the last two lines of Eq. (27), the plane wave $e^{i\bm{G}\cdot\bm{r}}$ was expanded in terms of spherical harmonics $Y_{lm}$ and spherical Bessel functions $j_{l}$. If only the spherically-symmetric terms corresponding to $l=m=0$ are kept, Eq. (27) simplifies to the spherically-symmetric atomic potential given in Eq. (15).

Appendix D Semilocal and non-local

parts of Hamann’s pseudopotentials

For the atom $\kappa$, the total non-local Vanderbilt-Kleinman-Bylander (NL-VKB) [56, 57, 58] pseudopotential corresponding to the angular momentum $l$ is generally written as

[TABLE]

where $V^{l}_{\kappa,\text{NL}}(r,r^{\prime})$ is the fully non-local part that includes only the beta-projectors $\beta_{\kappa,i}(r)$ and $\beta_{\kappa,j}(r)$ generated for the angular momentum $l$. These projectors are read directly from the pseudopotential file along with the normalization coefficients $D_{\kappa,ij}$. The local part $V_{\kappa,\text{loc}}(r)$ is extracted from a DFT calculation. With fully non-local pseudopotentials, the radial Schrödinger equation for the radial functions is

[TABLE]

While Eq. (29) can be solved numerically, in some cases it is more straightforward to convert the non-local pseudopotential to a semilocal form, if such a conversion is possible. The non-local Hamann’s ONCVPSP pseudopotentials [59, 60], used in the present work, allow for this conversion, as the pseudopotential files include the radial functions $\chi_{\kappa,l}(r)=rR_{\kappa,l}(r)$, evaluated at the same refernce energies used to construct the beta-projectors. With these radial functions, the semilocal components of the potential can be reconstructed as

[TABLE]

Alternatively, since the Hamann’s ONCVPSP code computes the SL parts internally, they can be printed directly in the pseudopotential file and read from there. In our tests, the SL parts calculated with Eq. (30) were found to be identical to those printed directly by the Hamann’s ONCVPSP code.