（常规）超导转变温度的第一性原理计算#

以 fcc-Pb 为例，使用 Quantum ESPRESSO 的 PHonon 模块计算电声耦合，并通过 McMillan–Allen–Dynes 公式估计常规超导转变温度。介绍结构优化、两次 SCF、声子与电声耦合计算、数据后处理和常见问题。

目的：在这里，我们根据DFPT理论和Eliashberg理论，计算常规超导体系中的电声耦合（electron-phonon coupling，EPC，有时也记为elph），并根据经验公式——McMillan公式估计超导转变温度。

计算软件：Quantum ESPRESSO（QE）中的PHonon模块，源码在QE/PHonon/PH路径

理论：BCS+Eliashberg理论

[1]Flores-Livas J A, Boeri L, Sanna A, et al. A perspective on conventional high-temperature superconductors at high pressure: Methods and materials[J]. Physics Reports, 2020, 856: 1-78.

[2]高淼, 孔鑫, 卢仲毅, 等. Li2C2 中电声耦合及超导电性的第一性原理计算研究[J]. 物理学报, 2015, 64(18): 186301.

Eliashberg 谱函数

http://www.w3.org/1998/Math/MathML" display="block"> α^{2} F (ω) = \frac{1}{2 π N (E_{F})} \sum_{q ν} \frac{γ_{q ν}}{ℏ ω_{q ν}} δ (ω - ω_{q ν})

声子线宽 linewidth

http://www.w3.org/1998/Math/MathML" display="block"> γ_{q ν} = \frac{2 π ω_{q ν}}{Ω_{B.Z}} \sum_{i, j} \int d^{3} k {| g_{k, q ν}^{i, j} |}^{2} δ (ε_{i, q} - E_{F}) δ (ε_{j, k + q} - E_{F})

电声耦合系数

http://www.w3.org/1998/Math/MathML" display="block"> λ = \sum_{q ν} \frac{γ_{q ν}}{π ℏ N (E_{F}) {ω_{q ν}}^{2}} = 2 \int_{0}^{\infty} \frac{α^{2} F (ω)}{ω} d ω

超导转变温度

http://www.w3.org/1998/Math/MathML" display="block"> T_{c} = \frac{f_{1} f_{2} ω_{\log}}{1.2} \exp (- \frac{1.04 (1 + λ)}{λ - μ^{*} (1 + 0.62 λ)})

对数声子频率

http://www.w3.org/1998/Math/MathML" display="block"> ω_{\log} = \exp (\frac{2}{λ} \int_{0}^{\infty} \frac{α^{2} F (ω)}{ω} \ln ω d ω)

计算基础：已经会用QE计算基本性质，包括结构优化、静态计算、（投影）能带结构/电子态密度，声子色散（声子谱）；如果不会，请参照教程：

https://yyyu200.github.io/DFTbook/blogs/2019/04/01/HandsOn/

接下来，我们以fcc-Pb为例，计算其超导转变温度。

赝势类型：ultrasoft pseudopotential（uspp）；泛函类型：PBE；

Pb.pbe-dn-rrkjus_psl.1.0.0.UPF

结果如下：

左图是声子色散，圆圈大小反映该mode对应的线宽linewidth；中间是声子态密度PHDOS；右边是超导谱函数和lambda曲线。根据MAD公式，估计Tc=4.3K。

前置内容#

需要强调的是，超导是金属在低温下的一种状态，所以只有金属才能计算超导。对于半导体/绝缘体或者半金属（semi-metal），没有超导电性一说。除非通过压力等方式使其相变，比如某结构常压下是bcc，在100GPa变成了fcc（通过焓压图+声子谱），材料由半导体变成了金属（能带计算），那就可以计算100GPa下fcc结构的超导转变温度。在这里，先验知识表明fcc-Pb是一个金属，其晶体结构（晶胞和原胞）和能带结构如下：

1. pw.x 结构优化#

pw.x做relax（或者vc-relax）：calculation=’vc-relax’

目的：得到0GPa在该赝势下最合理的晶格常数（不同赝势，不同泛函的晶格常数有区别），同时将结构内部原子的力释放（否则结构可能不稳定，有虚频）。

kmesh取18 18 18，是根据vaspkit的0.02*2π的密度生成，也可以取0.03

2. 第一次 SCF：生成致密 k 点数据#

pw.x做第一次scf计算

操作：将relax.out中优化完的结构提取，并写入scf0.in文件。设置calculation=’scf’，并加入参数la2F = .true.，这一步的kmesh取36 36 36（更密一些）

目的：设置la2F=true，会在./tmp目录下生成一个Pb.a2Fsave文件，这个文件是后面计算epc矩阵元所需要的。scf0.in文件如下：

&control
  calculation = 'scf',
  prefix = 'Pb',
  pseudo_dir = '/public/home/zpwu/QE/potential/us+pbe',
  outdir = './tmp',
  tprnfor = .true.,
  tstress = .true.,
  etot_conv_thr = 1.0d-7,
  forc_conv_thr = 1.0d-5,
/
&system
  ibrav = 0,
  nat = 1, ntyp = 1,
  occupations = 'smearing',
  smearing = 'mp',
  degauss = 0.02,
  ecutwfc = 80, ecutrho = 800,
  la2F = .true.,
/
&electrons
  electron_maxstep = 100,
  mixing_beta = 0.7,
  conv_thr = 1.0d-9,
/
&ions
  ion_dynamics = 'bfgs',
/
&cell
  cell_dofree = 'all',
  cell_dynamics = 'bfgs',
  press = 0,
  press_conv_thr = 0.01,
/
ATOMIC_SPECIES
Pb 207.2 Pb.pbe-dn-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (angstrom)
-2.4401838502  0.0000000000  2.4401838502
 0.0000000000  2.4401838502  2.4401838502
-2.4401838502  2.4401838502  0.0000000000
ATOMIC_POSITIONS (crystal)
Pb 0.0000000000 0.0000000000 0.0000000000
K_POINTS {automatic}
36 36 36 0 0 0

3. 第二次 SCF：生成声子计算波函数#

pw.x做第二次scf计算

操作：在将scf0.in文件cp为scf.in，并把la2F = true那一行注释掉，kmesh由36改成18

目的：重新生成波函数文件，用于下一步的ph.x计算。在这里注释掉la2F是避免覆盖掉原先致密k点生成的Pb.a2Fsave文件。

4. ph.x 计算电声耦合#

ph.x计算电声耦合

建立一个ph.in文件内容如下：

Electron-phonon coefficients
&inputph
  tr2_ph = 1.0d-15,
  prefix = 'Pb',
  fildvscf = 'Pbdv',
  outdir = './tmp',
  fildyn = 'Pb.dyn',
  electron_phonon = 'interpolated',
  el_ph_sigma = 0.005,
  el_ph_nsigma = 10,
  trans = .true.,
  ldisp = .true.,
  nq1 = 6, nq2 = 6, nq3 = 6
/

红框的参数就是控制计算EPC的，qmesh=666（nq1 nq2 nq3），对应kmesh=18 18 18。Kmesh和qmesh直接决定了ph.x的计算量。原则上qmesh要稍大一些，一般取kmesh的因子较为合适（比如一半或者三分之一）。

需要注意的是：这里的prefix参数和outdir要和pw.x计算时的in文件保持一致，以便于程序能顺利读取文件。

elph.sigma参数的物理意义：

高斯展宽近似

http://www.w3.org/1998/Math/MathML" display="block"> I = \int d k \int d q f (k, q) δ (ε_{k} - E_{F}) δ (ε_{k + q} - E_{F}) ≃ \frac{{Ω_{BZ}}^{2}}{N_{k} N_{q}} \sum_{k} \sum_{q} f (k, q) \frac{1}{\sqrt{2 π} σ} \exp (- \frac{{(ε_{k} - E_{F})}^{2}}{σ^{2}}) \times \frac{1}{\sqrt{2 π} σ} \exp (- \frac{{(ε_{k + q} - E_{F})}^{2}}{σ^{2}})

程序在计算delta函数时，通过近似成gauss函数来实现，因此引入了展宽broadening，当展宽足够小时，gauss函数近似为delta函数

其余参数的含义：https://www.quantum-espresso.org/Doc/INPUT_PH.html

运行ph.x开始计算EPC：mpirun -np Cores path/ph.x < ph.in > ph.out

mpirun -np Cores path/ph.x < ph.in > ph.out

这条命令写进任务提交脚本，Cores就是核数，path是ph.x所在路径，一般是QE的bin目录。如果调用的核数较多（32个），可以在ph.x后面写ph.x -nk 4 <ph.in >，加入nk参数来加速计算，nk后面的数字取不可约k点（scf.out文件）和核数的因子即可。

ph.x这一步是最耗时的（大约是pw.x的100~10000倍，取决于qmesh），等任务运行一两分钟后，可以查看ph.out文件的信息，在前面有这些内容：

Dynamical matrices for ( 6, 6, 6) uniform grid of q-points
(16 q-points):
 N      xq(1)         xq(2)         xq(3)
 1   0.000000000   0.000000000   0.000000000
 2  -0.117851130   0.117851130  -0.117851130
 3  -0.235702260   0.235702260  -0.235702260
 4   0.353553391  -0.353553391   0.353553391
 5   0.000000000   0.235702260   0.000000000
 6  -0.117851130   0.353553391  -0.117851130
 7   0.471404521  -0.235702260   0.471404521
 8   0.353553391  -0.117851130   0.353553391
 9   0.235702260   0.000000000   0.235702260
10   0.000000000   0.471404521   0.000000000
11   0.589255651  -0.117851130   0.589255651
12   0.471404521  -0.000000000   0.471404521
13   0.000000000  -0.707106781   0.000000000
14   0.471404521  -0.235702260   0.707106781
15   0.353553391  -0.117851130   0.589255651
16  -0.235702260  -0.707106781   0.000000000

6x6x6=216个q点，根据对称性，只有16个不等价的，于是ph.x只计算这16个q点，相应的坐标是在xyz基失下的坐标以2π/alpha为单位。

随后会开始进入第一个q点（0,0,0）的EPC计算（下图）。先是根据迭代计算该q点处每个mode的频率，由于对称性，Gamma点三个modes是等价的，因此只计算一个representation即可。ddv_scf2就是迭代的收敛值，当该值小于ph.in中的tr2_ph时，迭代结束；接着计算下一个representation.

Representation # 1 modes # 1 2 3

Self-consistent Calculation

Pert. # 1: Fermi energy shift (Ry) = -5.5145E-25   1.5673E-36
Pert. # 2: Fermi energy shift (Ry) =  4.4806E-25   2.6331E-36
Pert. # 3: Fermi energy shift (Ry) =  1.4648E-25  -9.4040E-38

iter # 1 total cpu time: 3.6 secs av.it.: 8.9
thresh= 1.000E-02 alpha_mix = 0.300 |ddv_scf|^2 = 5.253E-07

Gamma点计算顺利完成时，会输出该q点的频率（下图），以及相应的EPC-gamma line（GHz单位的），还有该mode对应的电声耦合系数lambda。由于Gamma点的声学支代表体系的平动，因此这里的gamma-linewidth和lambda为0。相应地会计算10个gaussian broadening下的lambda，对应ph.in里的sigma=0.005（单位是Ry）和nsigma=10。注意到这里的频率是负值-7cm-1，这是计算误差，可以忽略，后面在q2r.x中可以消除。但如果这里的虚频达到-20cm-1或者-50cm-1，则很可能虚频来源于光学支，意味着体系不稳定（畸变）。

Number of q in the star = 1
List of q in the star:
1 0.000000000 0.000000000 0.000000000

Diagonalizing the dynamical matrix
q = ( 0.000000000 0.000000000 0.000000000 )
freq (1) = -0.226702 [THz] = -7.561959 [cm-1]
freq (2) = -0.226702 [THz] = -7.561959 [cm-1]
freq (3) = -0.226702 [THz] = -7.561959 [cm-1]

electron-phonon interaction ...
Gaussian Broadening: 0.005 Ry, ngauss= 0
DOS = 4.203587 states/spin/Ry/Unit Cell at Ef= 11.915649 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz
Gaussian Broadening: 0.010 Ry, ngauss= 0
DOS = 3.595382 states/spin/Ry/Unit Cell at Ef= 11.920579 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz

由于一共有16个q点，因此这样的计算会重复16次。由于对称性，不同q点的计算量会不一样（Gamma点最快）。

如果计算顺利完成，则当前目录下会有16+1个Pb.dyn文件，一个是Pb.dyn0，里面就是16个不可约q点的坐标。打开Pb.dyn1，里面是动力学矩阵元信息以及频率信息。除此之外，还有一个elph_dir的文件夹，里面有一堆文件，我们关心elph.inp_lambda.N文件，一共有16个，对应每个q点的EPC，可以自行查看。比如第一个：

0.000000 0.000000 0.000000 10 3
-0.474854E-08 -0.474854E-08 -0.474854E-08
Gaussian Broadening: 0.005 Ry, ngauss= 0
DOS = 4.203587 states/spin/Ry/Unit Cell at Ef= 11.915649 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz
Gaussian Broadening: 0.010 Ry, ngauss= 0
DOS = 3.595382 states/spin/Ry/Unit Cell at Ef= 11.920579 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz
Gaussian Broadening: 0.015 Ry, ngauss= 0
DOS = 3.374842 states/spin/Ry/Unit Cell at Ef= 11.922760 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz
Gaussian Broadening: 0.020 Ry, ngauss= 0
DOS = 3.335246 states/spin/Ry/Unit Cell at Ef= 11.923377 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz
lambda( 2)= 0.0000  gamma= 0.00 GHz
lambda( 3)= 0.0000  gamma= 0.00 GHz
Gaussian Broadening: 0.025 Ry, ngauss= 0
DOS = 3.342013 states/spin/Ry/Unit Cell at Ef= 11.922985 eV
lambda( 1)= 0.0000  gamma= 0.00 GHz

第一行是q点坐标，10组展宽，3个modes（原胞1个原子，所以有3个modes）

第二行是频率信息，和dyn文件的频率是对应的，有个公式关系，好像对这个-0.4748E-8开个根号再乘以一个因子就到了dyn的频率，具体形式忘记了。这里是负值，表示虚频，不过Gamma点的小虚频不用管。

第三行开始就是10组展宽下的gamma linewidth和lambda值。

5. q2r.x 变换力常数#

运行q2r.x，对力常数做变换

需要一个q2r.in文件，内容如下：

&input
  zasr = 'simple',
  fildyn = 'Pb.dyn',
  flfrc = 'Pb.fc',
  la2F = .true.
/

这里zasr就是选取求和规则，这里选用simple

参数含义见：https://www.quantum-espresso.org/documentation/input-data-description/

和前面ph.x一样，脚本中加入mpirun …… path/bin/q2r.x < q2r.in > q2r.out，并提交任务。

mpirun …… path/bin/q2r.x < q2r.in > q2r.out

运行成功会有Pb.fc以及dyna2F文件，暂时用不到这两个文件，可以不管。

6. matdyn.x 计算声子色散与线宽#

运行matdyn.x根据力常数得到高对称点路径的声子色散关系及线宽

准备一个matdyn0.in文件如下：

&input
  asr = 'simple',
  flfrc = 'Pb.fc',
  flfrq = 'Pb.freq',
  flvec = 'Pb.modes',
  fleig = 'Pb.eig',
  la2F = .true.,
  q_in_band_form = .true.,
  q_in_cryst_coord = .true.,
  dos = .false.  ! for phonon band
/
8
0.0000000000 0.0000000000 0.0000000000 20 ! GAMMA
0.5000000000 0.0000000000 0.5000000000 20 ! X
0.6250000000 0.2500000000 0.6250000000 20 ! U
0.3750000000 0.3750000000 0.7500000000 20 ! K
0.0000000000 0.0000000000 0.0000000000 20 ! GAMMA
0.5000000000 0.5000000000 0.5000000000 20 ! L
0.5000000000 0.2500000000 0.7500000000 20 ! W
0.5000000000 0.0000000000 0.5000000000 20 ! X

下面的8个点就是倒格子的高对称点路径，在倒格式b1b2b3基失下的分数坐标，和能带路径一致即可。数字20表示每两个q点之间撒点20个，取50也可以。注意：这里的倒格子路径要和scf.in里面的结构对应。因为这里的scf.in是取的fcc结构原胞，其倒空间对称点就是图中这些。如果是取的fcc晶胞（四个原子）计算的，那倒空间对应的是simple cubic的路径，即G-X-M-R这些点。可查看：Setyawan W, Curtarolo S. High-throughput electronic band structure calculations: Challenges and tools[J]. Computational materials science, 2010, 49(2): 299-312.

参数含义见：https://www.quantum-espresso.org/documentation/input-data-description/

脚本提交：mpirun …… path/bin/matdyn.x < matdyn0.in > matdyn0.out

mpirun …… path/bin/matdyn.x < matdyn0.in > matdyn0.out

计算结束会输出10个elph.gamma.N文件，里面是每个展宽下高对称点路径上每个mode对应的gamma linewidth，后面我们会对这个文件的数据做后处理。

此外，还有名为Pb.freq和Pb.freq.gp的文件，里面是声子谱的信息，我们可以直接用.gp文件画图（下载到本地电脑，把后缀由gp改为dat或者txt，然后origin等软件画图即可）

还有一个Pb.modes的文件，里面是振动模式，可以借助工具进行可视化：

https://interactivephonon.materialscloud.io/

7. matdyn.x 计算谱函数与声子态密度#

运行matdyn.x，对qmesh网格插值得到电声作用谱函数和声子态密度

准备一个matdyn.in文件

&input
  asr = 'simple',
  flfrc = 'Pb.fc',
  la2F = .true.,
  dos = .true.,
  nk1 = 10, nk2 = 10, nk3 = 10,
  fldos = 'phonon.dos',  ! phonon dos
  ndos = 1000
/

这里nk1是对qmesh插值的精度，一般10~20足够，建议取qmesh的2~3倍（如果是二维材料，qmesh=661，那这里可以取12 12 1，z方向不用），太大了计算量大，而且精度提升不明显。

dos=true时，matdyn会计算声子态密度和电声作用谱函数；

ndos=1000，频率的步长，一般1000~2000足够。

注意，这里没有对flfreq等参数命令，因此程序会取默认值matdyn.freq。避免覆盖掉前一步matdyn.x生成的Pb.freq等文件。

脚本提交：mpirun …… path/bin/matdyn.x < matdyn.in > matdyn.out

mpirun …… path/bin/matdyn.x < matdyn.in > matdyn.out

最后会输出声子态密度文件phonon.dos，第一列是频率，单位是cm-1，第二列是total-phdos，第三列开始是每个原子对应的projected-phdos，按照scf.in里的原子顺序。下载到本地电脑，把dos后缀改成dat或者txt再用画图软件画。声子态密度的积分等于原子数N*3（无论是二维还是三维体系）。由于声子态密度做投影时的基矢是正交归一的（与投影电子dos不同），因此所有原子投影phdos之后等于total-phdos。

还有10个a2F.dos文件，对应10组展宽；a2F.dos文件第一列是频率，单位是Ry，第二列是谱函数a2F(w)，无量纲，第三列开始是每条声子模式对应的谱函数，这里原胞1个原子，对应三条模，因此后面有三列

# Eliashberg function a2F (per both spin)
# frequencies in Rydberg
# DOS normalized to E in Rydberg: a2F_total, a2F(mode)
0.412614E-06 0.174960E-09 0.124253E-09 0.490829E-10 0.162485E-11
0.123785E-05 0.472395E-08 0.354482E-08 0.132525E-08 0.438737E-10
0.206309E-05 0.218702E-07 0.155316E-07 0.613545E-08 0.203122E-09
0.288834E-05 0.600117E-07 0.426187E-07 0.168357E-07 0.557369E-09
0.371358E-05 0.127547E-06 0.905802E-07 0.357820E-07 0.118462E-08
0.453882E-05 0.232874E-06 0.165380E-06 0.653304E-07 0.216286E-08
0.536406E-05 0.384390E-06 0.272983E-06 0.107837E-06 0.357011E-08
0.618930E-05 0.590494E-06 0.419353E-06 0.165657E-06 0.548435E-08
0.701454E-05 0.859585E-06 0.610453E-06 0.241148E-06 0.798360E-08
0.783978E-05 0.120006E-05 0.852249E-06 0.336665E-06 0.111458E-07
0.866502E-05 0.162032E-05 0.115070E-05 0.454564E-06 0.150491E-07

一个lambda文件，里面是10组展宽下的电声耦合系数和对数频率

Electron-phonon coupling constant, lambda

Broadening   lambda  dos(Ef)  omega_ln [K]
0.0050       1.2839   4.2036    75.0949
0.0100       1.0139   3.5954    72.6999
0.0150       0.9121   3.3748    71.6661
0.0200       0.8883   3.3352    71.0829
0.0250       0.8685   3.3420    70.7954
0.0300       0.8894   3.3557    70.6804
0.0350       0.8919   3.3689    70.6385
0.0400       0.8931   3.3809    70.6276
0.0450       0.8927   3.3918    70.6372
0.0500       0.8904   3.4016    70.6671

最后需要根据lambda文件，选定一组合适的展宽，在这里选的是第四组，broadening=0.02Ry。查看a2F.dos4，最后一行同样会输出由谱函数所计算得到的lambda，

0.824828E-03 0.000000E+00 0.000000E+00 0.000000E+00
lambda = 0.882761851404398    Delta = 8.252406453807191E-007

二者的差别源于计算精度。这个展宽该如何选取，我们最后会单独讨论

8. 数据后处理与画图#

数据后处理+画图

（1）要得到Tc，需要根据a2F.dos4文件，通过一个小脚本alpha_lambda.py，计算得到epc系数和对数频率，并根据McMillan-Allen-Dynes公式得到Tc，储存在modes.txt文件前几行。

（2）同时也会输出单位换算后的谱函数与lambda曲线数据a2f_lambda.txt

（3）对Pb.freq和elph.lambda.4的数据进行处理，得到方便画图的数据格式。原理：这两个文件的q点都是一一对应的，画图的时候把Pb.freq化成散点图，点的大小以elph.lambda.4中每个phonon mode的linewidth数据作为权重。由于两个文件并非Pb.freq.gp那种可直接画图的格式，因此需要另外的脚本后处理。这里采用lambda_q.py脚本。最后会输出一个final.txt文件，第一列是序号，与matdyn0.in中的q-path对应。第二列是频率，cm-1单位；第三列是对应线宽linewidth；第四列是对应的lambda_qv，公式在第一页。

将数据最终呈现如下：

附：通常高对称点之间对应在x轴的长度是不同的，在这里我们一律取成一样长了，这不影响结论。

附录：常见问题#

1. lambda 展宽的选取#

（1）关于lambda展宽的选取，一般要求broadening不能太大，经验上在0.02~0.06Ry量级是合适的。主要是因为qmesh没法取得足够密，比如取到12x12x12开始计算量就很大了。为了解决这个问题，后来发展出了一套基于Wannier函数的方法，是EPW软件，属于QE的一个模块。但EPW其实也没法很好地解决这一问题。

可以详细看：

[1]Wierzbowska M, de Gironcoli S, Giannozzi P. Origins of low-and high-pressure discontinuities of Tc in niobium[J]. arXiv preprint cond-mat/0504077, 2005.

[2]Shipley A M, Hutcheon M J, Johnson M S, et al. Stability and superconductivity of lanthanum and yttrium decahydrides[J]. Physical Review B, 2020, 101(22): 224511.

[3]Gao M, Li Q Z, Yan X W, et al. Prediction of phonon-mediated superconductivity in borophene[J]. Physical Review B, 2017, 95(2): 024505.

2. 库伦屏蔽常数的选取#

（2）关于库伦屏蔽常数 $http://www.w3.org/1998/Math/MathML" display="inline"> μ^{*}$ 的选取，一般取经验值0.1~0.13左右。原则上和体系有关，也有一些从第一性原理出发估计该值的。

Morel P, Anderson P W. Calculation of the superconducting state parameters with retarded electron-phonon interaction[J]. Physical Review, 1962, 125(4): 1263.

Bauer J, Han J E, Gunnarsson O. Retardation effects and the Coulomb pseudopotential in the theory of superconductivity[J]. Physical Review B, 2013, 87(5): 054507.

3. 虚频与超导计算#

（3）如果声子谱有虚频，可以计算超导吗？

分情况讨论。如果是Gamma点附近的虚频，且可以验证是计算精度导致的虚频，并且该频率对epc贡献很弱（通过分析linewidth和对应的lambda），那在计算超导时可以直接把虚频部分截掉，谱函数积分时只对正频率那些数据计算即可。

如果线宽数据显示虚频或者接近0频率的声子对epc贡献很大，比如谱函数在0频率处有异常高的峰，导致lambda曲线在接近0的低频处有很高的平台，那说明虚频对epc的影响不可忽略。

4. qmesh 收敛性测试#

（4）需要对qmesh做收敛性测试吗？

需要，一般需要计算不同qmesh下的Tc，原则上qmesh要大一点。在算力允许的前提下选能达到Tc/lambda收敛的qmesh。个人经验，在kmesh取到2πx0.025的密度下，q-mesh取kmesh的一半甚至1/3是不错的选择。

5. Tc 计算精度#

（5）QE/PHonon能算准Tc吗？

不好说，第一性原理计算想彻底算准某个性质还是具有挑战性的，但至少通常可以给出不错的结果。

6. 多任务并行计算 q 点#

（6）ph.x模块计算16个不可约q点太慢了，是否有更高效的办法？

有，可以通过多个任务同时进行，同时计算不同q点。因为不同q点之间的计算是独立的，因此我们可以在a任务中计算1-4个q点，b任务计算5-8个q点，将16个q点分成4个任务同时进行。（也可以分成8个，取决于你可调用的计算资源）

通过在ph.in中设置start_q = 1, last_q = 4，可以控制Ph.x只计算第1-4个q点。

记得要新建不同的文件夹，在里面计算。可以从scf0.in开始重新算。

计算完成后，把所有dyn文件复制到同一个目录下，并且把abcd四个路径下elph_dir文件夹中的文件全部复制到同一个elph_dir中，要保证16+1个dyn文件和elph_dir文件夹是在同一路径的。然后开始q2r.x的计算。

7. ph.x 不收敛#

（7）ph.x计算不收敛怎么办？

首先确认收敛情况，是直接发散，还是在某个值附近振荡。可以查看ph.out文件确认

Representation # 1 modes # 1 2 3

Self-consistent Calculation

Pert. # 1: Fermi energy shift (Ry) = -5.5145E-25   1.5673E-36
Pert. # 2: Fermi energy shift (Ry) =  4.4806E-25   2.6331E-36
Pert. # 3: Fermi energy shift (Ry) =  1.4648E-25  -9.4040E-38

iter # 1 total cpu time: 3.6 secs av.it.: 8.9
thresh= 1.000E-02 alpha_mix = 0.300 |ddv_scf|^2 = 5.253E-07

如果|ddv_scf|2的值一直在1E-10或者1E-12附近，离预设值比如1E-14有一些距离，那么：

a.在ph.in里面设置alpha_mix=0.3来改善收敛情况（默认alpha_mix=0.7），也可以尝试其它值。

b.或者增加ph.x的最大步数，通过修改PH源码参数实现，自行搜索ph.x的niter_ph参数

c.网上搜一搜

8. q 点坐标的基矢#

（8）不同文件中q点坐标基矢区别

Xxx.dyn0和一系列dyn里面的q点坐标是xyz基矢下的坐标，在dynN文件里有如下说明：

Dynamical Matrix in cartesian axes
q = ( 0.000000000 -0.235702527 0.333330313 )

1 1
0.46048770 0.00000000 -0.00000000 0.00000000 0.00000000 -0.00000000
0.00000000 0.00000000 0.46158748 0.00000000 -0.06313408 0.00004308
0.00000000 0.00000000 -0.06313408 -0.00004308 0.71604087 0.00000000

ph.out里面的q点坐标也是xyz基矢下的分数坐标

在第一次matdyn.x运行后输出的xxx.freq文件和xxx.modes文件里面，q点坐标是在xyz基矢下的坐标

Matdyn0.in输入的声子路径是在b1b2b3基失下的

例如：

某结构的matdyn0.in里面路径如下

0.000000  0.000000 0.000000 40 !$\Gamma$
0.500000  0.000000 0.000000 40 !L
0.500000 -0.500000 0.500000 40 !T
0.251363  0.251363 0.251363 40 !P$_2$
0.000000  0.000000 0.000000 40 !$\Gamma$
0.500000 -0.500000 0.000000  1 !F

这些q点是在b1b2b3基矢下的分数坐标，我们讨论最后一个点F，即(0.5,-0.5,0)

再查看ph.out或者scf.out文件前面几百行输出的倒格矢信息如下

crystal axes: (cart. coord. in units of alat)
a(1) = ( 0.816496  0.471404 0.333336 )
a(2) = (-0.816496  0.471404 0.333336 )
a(3) = ( 0.000000 -0.942808 0.333336 )

reciprocal axes: (cart. coord. in units 2 pi/alat)
b(1) = ( 0.612373  0.353554 0.999991 )
b(2) = (-0.612373  0.353554 0.999991 )
b(3) = ( 0.000000 -0.707108 0.999991 )

由此可以推得F点在xyz基矢下的坐标为

0.5b1 – 0.5b2 +0 = (0.612373, 0 0 )

此时我们查看xxx.freq文件最后一个q点

942.7019 958.3146 1001.3767
0.612373 0.000000 0.000000
-21.9767 133.1106 197.3958 236.8210 275.0763 284.4752
296.4678 325.5081 343.5939 368.6282 381.1751 382.9941
397.1765 407.5584 410.6235 431.7804 447.0141 447.8363
455.3328 468.0262 472.9859 480.5742 491.3800 497.6897
504.2299 509.4128 516.4381 529.9273 540.3243 547.9220
550.4473 551.1212 556.2235 561.6376 568.9302 574.3356
578.2211 590.0323 597.8327 602.8434 602.8899 618.1930
621.7575 631.3807 637.7752 640.2146 642.8590 664.5740

与我们的推断相符

需要说明的是pw.x输出的out文件前面几十行关于k点坐标的信息中，k点坐标也是xyz基失下的，如下

number of k points= 6  Methfessel-Paxton smearing, width (Ry)=0.0200
cart. coord. in units 2pi/alat
k(1) = ( 0.0000000  0.0000000 0.0000000), wk = 0.0740741
k(2) = ( 0.0000000 -0.2357025 0.3333303), wk = 0.4444444
k(3) = (-0.2041244 -0.1178513 0.6666606), wk = 0.4444444
k(4) = (-0.2041244  0.3535538 0.0000000), wk = 0.4444444
k(5) = ( 0.0000000  0.0000000 0.9999909), wk = 0.1481481
k(6) = ( 0.0000000  0.4714051 0.3333303), wk = 0.4444444

9. 虚频排查#

（9）关于虚频

算得的dyn文件里面有虚频（Gamma点之外），那声子谱也必然有虚频。因为声子谱是根据这些dyn的数据插值得来。需要通过提高计算精度排查虚频是否“为真”。

（1）提高ph.in里面的tr2_ph=1d-15~19

（2）确保算声子的k点密度不低于2πx0.03甚至0.025（个人经验规律，原则上要自己测试频率随k点变化）

（3）确保截断能和电荷密度足够，至少要保证scf的收敛性。USPP和PAW一般70Ry以上足够，ONCV一般90~100Ry。对于一些稀土f元素和部分d元素，这个值可能更大

（4）scf.in里面的deguass取大一些，有利于消除dyn里面的虚频。Smearing如果取的guassian或者mp，deguass可以取到0.02甚至0.04Ry。但太大对应的结果就不合理了，因为这个值控制电子在Ef附近的占据情况，太高的deguass对应着“电子温度”很高。一般只有高温AIMD才会用到高deguass值（~0.05Ry）

如果dyn里面没有虚频，但是matdyn.x之后得到的声子谱有虚频

更换q2r.in和matdyn.in的求和规则，simple改成crystal（crystal的计算量大很多，甚至对某些原子多的体系~20atoms可能会内存不够报错）

可能因为ph.in里面q-mesh采样不够导致。检验办法是单独计算虚频处q点的频率

对应的ph.in文件如下

Electron-phonon coefficients
&inputph
  tr2_ph = 1.0d-14,
  prefix = 'pwscf',
  fildvscf = 'YOUR_DVSCF',
  outdir = './tmp',
  fildyn = 'YOUR_DYN',
  electron_phonon = 'interpolated',
  el_ph_sigma = 0.005,
  el_ph_nsigma = 10,
  alpha_mix = 0.3,
  !trans = .true.,
  !ldisp = .true.,
  !start_q = 1,
  !last_q = 1,
  !nq1 = 3, nq2 = 3, nq3 = 3
  search_sym = .false.,
  nmix_ph = 8
/
0.612373 0.000000 0.000000

请将 YOUR_DVSCF、YOUR_DYN 替换为实际文件名；search_sym = .false.。

注释了ldisp参数和nq参数，并在最后/后面一行写了相应q点坐标（这个坐标是xyz基矢下的，可以从xxx.freq文件里面找到）

如果算得的q点有虚频，那重复上一步的虚频检验方法；如果没有，说明需要取更密的q-mesh来得到无虚频的声子谱。

注：

Dyn文件里面的本征矢量是在a1a2a3基失下的

Matdyn输出的modes后缀文件的本征矢量是在xyz基失下的

在vesta标定振动模式时，取dyn方便