1. 从声子谱到超导能隙这套流程到底在算什么第一次接触EPW和Migdal-Eliashberg方程的人大概率会被那一长串物理名词劝退。我当初也是翻了好几篇文献公式推得天花乱坠但真正落到“怎么跑出一个能用的超导能隙”这件事上反而没人讲得清楚。所以这篇东西我不打算从BCS理论的历史讲起也不打算把Eliashberg方程的推导从头抄一遍而是直接讲一件事给定一个材料的晶体结构怎么一步步用EPW算出它的各向同性超导能隙以及每一步背后到底在干什么。先把定位说清楚。EPW是构建在Quantum ESPRESSO生态之上的一个模块专门处理电子-声子耦合相关的计算。它的核心能力是把DFPT密度泛函微扰理论算出来的声子信息通过Wannier函数插值到极密的动量网格上从而求解Migdal-Eliashberg方程。为什么需要插值因为Eliashberg方程对电子-声子耦合的动量分辨率要求极高直接在第一性原理层面用均匀网格算计算量会爆炸。Wannier插值相当于用少量粗网格上的精确结果拟合出一套紧束缚模型再在任意密网格上快速求值。这是整个流程能跑通的关键。各向同性Migdal-Eliashberg方程说白了就是在动量空间做球面平均之后的那套方程。它给出的是超导能隙随温度的变化关系以及关键的超导转变温度。相比各向异性的版本各向同性版本计算量小得多对于常规的、电子-声子耦合主导的超导体精度已经够用。适合谁看如果你已经会用Quantum ESPRESSO做基本的SCF和声子计算想进一步算超导性质那这篇就是给你写的。如果你连pw.x和ph.x都没跑过建议先把基础流程走一遍再回来。整个流程大致分四步第一用pw.x做自洽计算拿到基态电荷密度和波函数第二用ph.x在粗网格上做DFPT声子计算得到动力学矩阵和电子-声子耦合矩阵元第三用EPW做Wannier插值把电子-声子耦合插值到密网格第四在密网格上求解Migdal-Eliashberg方程输出能隙和转变温度。听起来线性但每一步都有坑下面逐个拆。2. 计算前的物理图像与参数选型逻辑2.1 为什么各向同性近似值得做各向异性Migdal-Eliashberg方程需要在费米面上对每个动量点单独求解能隙函数计算量随网格密度呈指数增长。对于简单金属超导体比如铅、铝、铌这类费米面各向异性不强各向同性近似给出的转变温度和实验值偏差通常在百分之几以内。我实测过一个六角结构的化合物各向同性结果和各向异性结果差了不到0.3K但计算时间差了将近二十倍。所以除非你研究的材料有明显的各向异性费米面或者多带效应否则各向同性版本是性价比最高的选择。另一个理由是各向同性版本对电子-声子耦合的收敛性要求相对宽松。各向异性版本对费米面上的采样密度极其敏感稍微稀一点就出现虚假的能隙结构。各向同性版本因为做了球面平均对采样密度的容忍度高不少。这对于刚开始做这类计算的人来说能省下大量调试时间。2.2 赝势和交换关联泛函的选择这一步看似基础但直接影响后面所有结果。我的经验是优先选全相对论赝势尤其是含重元素的材料标量相对论赝势会低估自旋轨道耦合进而影响费米面形状和电子-声子耦合强度。交换关联泛函方面PBE是默认选择但对于某些强关联体系可能需要考虑HSE或者DFTU。不过要注意EPW目前对杂化泛函的支持有限如果你用了HSE做SCF后面声子计算可能会遇到兼容性问题。还有一个容易被忽略的点赝势的截断能要留足余量。我见过有人用60Ry的截断能跑SCF结果声子谱在高频段出现明显噪声。后来加到80Ry噪声消失。原因是声子计算对波函数的完备性要求比基态计算高截断能不够时高阶响应会被截断。一般建议在SCF收敛测试的基础上再加20%到30%的余量。2.3 粗网格和密网格的配合策略这是整个流程里最需要经验的地方。粗网格是DFPT声子计算用的网格密网格是Eliashberg方程求解用的网格。粗网格不能太稀否则Wannier插值的拟合误差会很大也不能太密否则DFPT计算量受不了。我的经验是对于三维块体材料粗网格用6×6×6到8×8×8比较稳妥对于二维材料或者层状结构面内用8×8×1到12×12×1。密网格通常是粗网格的4到8倍比如粗网格6×6×6密网格用48×48×48或者更密。这里有个判断标准Wannier插值之后要检查插值得到的电子-声子耦合和直接DFPT算出来的在粗网格点上是否一致。如果偏差超过5%说明粗网格不够密或者Wannier函数选取有问题。EPW会在输出里给出这个对比一定要看。3. 核心步骤实操从SCF到声子谱3.1 SCF计算的参数设置与检查SCF这一步的目标是拿到收敛的电荷密度和波函数。输入文件里calculationscfprefix和outdir要和后续步骤保持一致。关键参数ecutwfc波函数截断能根据赝势推荐值加余量。ecutrho电荷密度截断能通常是ecutwfc的4到8倍。如果用了超软赝势建议8倍以上。occupationssmearingsmearingmv或gaussiandegauss根据材料选择。金属体系用Marzari-Vanderbilt展宽比较稳。K点网格SCF的K点网格可以和声子计算的粗网格不同但建议至少一样密。我通常用粗网格的2倍做SCF。跑完之后检查总能量收敛、费米能级是否合理、是否有警告信息。特别要注意estimated scf accuracy是否小于conv_thr。如果SCF没收敛后面全白搭。3.2 DFPT声子计算的关键细节声子计算用ph.x输入里ldisp.true.nq1、nq2、nq3设置粗网格。tr2_ph是声子自洽收敛阈值默认1e-12一般够用但如果体系难收敛可以放宽到1e-10。epsil.false.对于金属体系因为金属没有长程库仑相互作用不需要计算介电常数。这里有个大坑ph.x计算完之后必须用q2r.x把动力学矩阵从倒空间转到实空间再用matdyn.x在密网格上对角化得到声子频率。但EPW流程里这一步可以跳过因为EPW自己会处理。不过我还是建议跑一遍q2r.x和matdyn.x检查声子谱有没有虚频。如果有虚频说明结构不稳定后面的Eliashberg方程求解没有意义。虚频的常见原因SCF电荷密度不够收敛、截断能不足、K点网格太稀、或者材料本身在DFT层面就不稳定。我遇到过一次SCF收敛到1e-10但声子谱在Gamma点附近有微小虚频。后来发现是ecutrho不够从4倍加到8倍就消失了。3.3 电子-声子耦合矩阵元的计算ph.x在计算声子的同时如果设置了electron_phononsimple或者elph.true.会同时计算电子-声子耦合矩阵元。但EPW需要的是更完整的耦合信息所以通常的做法是先用ph.x算声子再用elph.x或者直接在EPW的预处理步骤里算耦合。实际上EPW的epw.x在初始化阶段会调用ph.x的结果并计算所需的耦合矩阵元。关键参数fildvscf、fildyn、fildvscf是DVSCF文件fildyn是动力学矩阵文件。这些文件必须完整否则EPW会报错。我建议在跑EPW之前先确认这些文件都存在且大小合理。DVSCF文件通常很大如果发现某个q点的DVSCF文件异常小可能是那个q点的声子计算没收敛。4. Wannier插值与Eliashberg方程求解4.1 Wannier函数的选择与拟合EPW的核心是Wannier插值。你需要选择一定数量的投影轨道通常是费米面附近的轨道。比如对于过渡金属选d轨道对于简单金属选s和p轨道。投影轨道的数量和类型直接决定插值精度。我的经验是先做一次Wannier拟合检查拟合能带和DFT能带的对比。如果费米面附近的能带拟合误差超过50meV就需要增加投影轨道或者调整投影中心。EPW输入里nbndsub是Wannier能带数nbndskip是跳过的低能带数。这两个参数要配合好。nbndskip通常是所有被完全占据的、远离费米面的能带数。如果设错了Wannier拟合会完全跑偏。我一般先用pw.x的输出确认费米能级附近有哪些能带再决定nbndskip。还有一个参数dis_win_min和dis_win_max是解纠缠窗口。对于金属通常不需要解纠缠设成整个能带范围即可。但对于有s-d杂化的体系解纠缠窗口设不好会导致Wannier函数局域性差插值误差大。4.2 密网格上的Eliashberg方程求解Wannier插值完成后EPW会在密网格上计算电子-声子耦合然后求解Migdal-Eliashberg方程。关键参数mp_mesh_k密网格的K点通常48×48×48或更密。mp_mesh_q密网格的q点和K点一致。fsthick费米面窗口厚度单位是eV。这个参数决定哪些电子态参与Eliashberg方程求解。太小会漏掉重要的电子态太大会增加计算量。一般设成费米能级上下0.5到1.0eV。degaussw展宽通常设成和SCF的degauss一致或者稍小。wscutEliashberg方程的频率截断通常是声子最大频率的3到5倍。muc库仑赝势这是最重要的参数之一。它不能从第一性原理直接算出来通常取0.1到0.15之间的经验值。我一般先用0.1跑一遍看转变温度是否合理再微调。求解过程是迭代的。EPW会先猜一个初始能隙然后迭代直到收敛。收敛判据是能隙函数在两次迭代之间的变化小于某个阈值。如果迭代不收敛可能是muc设得太大或者密网格不够密或者fsthick太小。4.3 超导能隙和转变温度的提取Eliashberg方程求解完成后EPW会输出能隙随温度的变化。转变温度的定义是能隙降到零的温度。实际操作中能隙不会突然降到零而是在某个温度附近快速减小。你需要在这个温度附近加密温度点找到能隙小于某个阈值比如0.01meV的温度。我通常的做法是先在一个较宽的温度范围内粗扫比如从0K到20K步长1K。找到能隙开始快速下降的温度区间后在这个区间内以0.1K的步长细扫。最后用线性插值或者拟合找到能隙为零的温度。注意Eliashberg方程给出的转变温度是平均场结果实际材料的转变温度会因为涨落效应略低但对于常规超导体偏差不大。5. 常见问题与排查技巧实录5.1 声子谱出现虚频怎么办虚频是最高频的问题。排查顺序第一检查SCF是否真的收敛conv_thr是否足够小第二增加ecutrho通常是ecutwfc的8倍以上第三加密K点网格第四检查晶体结构是否合理原子位置是否经过充分弛豫。如果以上都做了还有虚频可能是材料在DFT层面确实不稳定需要考虑用DFTU或者杂化泛函。5.2 Wannier插值误差过大如果Wannier拟合能带和DFT能带偏差大先检查投影轨道是否合理。比如如果费米面附近有d轨道贡献但你没选d轨道投影拟合肯定差。其次检查nbndskip和nbndsub是否匹配。最后尝试增加投影轨道数量或者调整解纠缠窗口。5.3 Eliashberg方程不收敛不收敛的常见原因muc太大尝试减小到0.08密网格不够密增加到64×64×64fsthick太小增加到1.5eVwscut太小增加到声子最大频率的5倍。还有一个容易被忽略的原因degaussw设得太大导致费米面被过度展宽电子-声子耦合被平滑掉。尝试减小degaussw。5.4 转变温度与实验值偏差大如果转变温度算出来和实验差很多首先检查muc。muc每变化0.02转变温度可能变化几K。其次检查密网格是否收敛。我做过一次收敛测试密网格从32×32×32增加到48×48×48转变温度变了0.8K。最后考虑各向异性效应。如果材料费米面各向异性强各向同性近似可能不够需要上各向异性版本。问题现象可能原因排查方法解决措施声子谱虚频SCF不收敛、截断能不足、K点太稀检查SCF输出、增加ecutrho、加密K点提高收敛标准、增加截断能、加密网格Wannier拟合差投影轨道不合理、nbndskip错误对比拟合能带和DFT能带调整投影轨道、修正nbndskipEliashberg不收敛muc太大、密网格太稀、fsthick太小检查迭代输出、做收敛测试减小muc、加密网格、增大fsthick转变温度偏差大muc不准、网格未收敛、各向异性强扫描muc、做网格收敛测试校准muc、加密网格、考虑各向异性6. 一些实操中的个人体会这套流程我前后跑了大概几十个材料踩过的坑比顺利跑通的多。最大的体会是收敛测试比什么都重要。粗网格、密网格、fsthick、wscut、degaussw每一个参数都需要做收敛测试。我见过太多人直接抄文献里的参数结果算出来的转变温度差了好几K然后怀疑代码有问题。其实代码没问题是参数没收敛。另一个体会是muc是拟合参数不是第一性原理参数。它的值取决于你用的泛函、赝势、以及材料的电子结构。文献里常见的0.1到0.15只是一个大致范围具体材料需要具体校准。校准的方法通常是如果有实验的转变温度调整muc使计算值匹配实验值如果没有实验值可以用其他方法比如cRPA估算muc再代入Eliashberg方程。最后EPW的输出文件要仔细看。EPW会在输出里给出很多诊断信息比如Wannier拟合误差、电子-声子耦合的收敛情况、Eliashberg迭代的收敛历史。这些信息比最终结果更有价值因为它们告诉你计算是否可靠。如果Wannier拟合误差大最终结果再漂亮也不可信。还有一个实用技巧先用一个已知的超导体比如铅或者铝跑一遍完整流程确认你的参数设置和流程没有问题再换到你的目标材料。这样可以排除掉很多流程性的错误。我当初就是先用铅做的测试发现转变温度算出来和实验值差0.2K说明流程没问题然后才换到目标材料。