电声耦合常数计算再谈2017

1) 电声耦合常数计算的基本流程

2) 氢化物超导机制探讨

3) 电声耦合常数计算公式背后的物理意义

4) 常见问题收集

mage              PC                Routine            Line        Source

lambda.x           00000000004A5BCD  Unknown               Unknown  Unknown

lambda太大了，输出都超格式了！

http://blog.sciencenet.cn/blog-567091-726265.html

scf fit 36 36 24
scf 18 18 12
q  6 6 4

A:::

1、在lambda.x的input文件的第一行，需要输入3个参量，分别是emax，degauss，以及smearing method。请问这三个参数有什么含义，取法有什么讲究，对于计算的结果有什么影响？

2、在lambda.x的output文件里面，分别列出了degauss从0.005到0.05的lambda及Tc的计算结果，不同的degauss对应的结果不同，那么，到底哪一个结果才是我们想要的最终结果呢

Q:::

emax是比最高声子能量还高一点的一个值，取多少自己定。声子能量，直接看最高声子频率就是，要取以THz为单位的。
degauss和smearing method是在布里渊区内取点的时候用的展宽方式和大小，但具体该怎么取我还不懂

，最后一行取0.1到0.15之间，是代表电子之间的库仑排斥作用

1) PWSCF如何有效建立晶体模型及查看建立模型的正确性：PWSCF建模与可视化---为了忘却的纪念

http://blog.sciencenet.cn/blog-567091-725337.html

http://web.mit.edu/xcrysden_v1.4.1/www/XCRYSDEN/doc/pwscf.html

ibrav = 0

CELL_PARAMETERS{ alat | bohr | angstrom }

v1(1)  v1(2)  v1(3)

v2(1)  v2(2)  v2(3)

v3(1)  v3(2)  v3(3)

如果直接利用该坐标系建立QE结构模型，则：
1） QE优化非常耗时；
2）优化得到的结果可能很诡异。
因此，正确的做法是：
将VASP优化后的结果，导入MS中施加对称性后，再通过vesta转化为VASP格式的坐标文件。此时，即可利用重新调整后的坐标系及原子坐标，进行QE建模。

&system
ibrav=0, celldm(1) = 1.8897268777743552,   ####直接将cell parameters中的A单位转化为Bohr
nat=14,
ntyp=2,

ATOMIC_POSITIONS (crystal)          (VASP中的原子坐标直接引用)
H    0.000000000         0.500000000         0.750000000
H    0.000000000         0.500000000         0.250000000
H     0.750000000         0.000000000         0.500000000
H     0.250000000         0.000000000         0.500000000
H     0.500000000         0.750000000         0.000000000
H     0.500000000         0.250000000         0.000000000
H     0.500000000         0.000000000         0.250000000
H     0.500000000         0.000000000         0.750000000
H     0.000000000         0.750000000         0.500000000
H     0.000000000         0.250000000         0.500000000
H     0.750000000         0.500000000         0.000000000
H     0.250000000         0.500000000         0.000000000
X    0.500000000         0.500000000         0.500000000
X     0.000000000         0.000000000         0.000000000
CELL_PARAMETERS (VASP中的坐标矢量直接引用）
3.2657001019         0.0000000000         0.0000000000
0.0000000000         3.2657001019         0.0000000000
0.0000000000         0.0000000000         3.2657001019

ecutwfc=60.0, ecutrho = 1000.0
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.05
/

2）如何获取超导温度

成功运行例子文件 example03  ：电声耦合常数计算参数设置  elph=.true.,      5.0.2版改为 electron_phonon='interpolated',

get-Tc的脚本

cat ../getlambda.sh

#!/bin/sh

cat > lambda.in << EOF

36.635825  0.01  1    ! emax (something more than highest phonon mode in THz), degauss, smearing method

16          ! Number of q-points for which EPC is calculated,

0.0000000   0.0000000   0.0000000  0.0092593

0.1178511   0.1178511  -0.1178511  0.0740741

0.2357023   0.2357023  -0.2357023  0.0740741

-0.3535534  -0.3535534   0.3535534  0.0370370

0.0000000   0.2357023   0.0000000  0.0555556

0.1178511   0.3535534  -0.1178511  0.2222222

-0.4714045  -0.2357023   0.4714045  0.2222222

-0.3535534  -0.1178511   0.3535534  0.2222222

-0.2357023   0.0000000   0.2357023  0.1111111

0.0000000   0.4714045   0.0000000  0.0555556

-0.5892557  -0.1178511   0.5892557  0.2222222

-0.4714045   0.0000000   0.4714045  0.1111111

0.0000000  -0.7071068   0.0000000  0.0277778

-0.7071068  -0.2357023   0.4714045  0.2222222

-0.5892557  -0.1178511   0.3535534  0.2222222

0.0000000  -0.7071068  -0.2357023  0.1111111

elph. 0.000000. 0.000000. 0.000000   ! elph output file names,

elph. 0.117851. 0.117851.-0.117851

elph. 0.235702. 0.235702.-0.235702   ! in the same order as the q-points before

elph.-0.353553.-0.353553. 0.353553

elph. 0.000000. 0.235702. 0.000000

elph. 0.117851. 0.353553.-0.117851

elph.-0.471405.-0.235702. 0.471405

elph.-0.353553.-0.117851. 0.353553

elph.-0.235702. 0.000000. 0.235702

elph. 0.000000. 0.471405. 0.000000

elph.-0.589256.-0.117851. 0.589256

elph.-0.471405. 0.000000. 0.471405

elph. 0.000000.-0.707107. 0.000000

elph.-0.707107.-0.235702. 0.471405

elph.-0.589256.-0.117851. 0.353553

elph. 0.000000.-0.707107.-0.235702

0.10                     ! \mu the Coloumb coefficient in the modified

! Allen-Dynes formula for T_c (via \omega_log)

EOF

/home/ah736/espresso/src/espresso-4.3.1/bin/lambda.x  < lambda.in > lambda.out

cat lambda.out

.........

lambda        omega_log          T_c

0.27157       487.422              0.079

0.23829       478.847              0.012

0.24074       475.229              0.014

0.24202       472.958              0.015

0.24165       470.570              0.015

0.24030       469.347              0.013

0.23834       468.010              0.012

0.23597       467.256              0.010

0.23352       466.806              0.008

0.23126       466.296              0.007

3） 从PWscf声子谱计算中获取零点振动能数据

http://blog.sciencenet.cn/blog-567091-727689.html

ZPE是物质的固有属性，不随计算方法而改变，即选取不同的参数如赝势种类、截断能大小等，只要精度足够，不影响计算得到的零点能数值。

4） PWSCF如何进行结构优化，各参数如何设置：pw.x    --run_optim脚本

5）PWSCF-计算声子谱

http://blog.sciencenet.cn/blog-567091-758948.html

number of k points=    20
grep "number"

6） 计算电声耦合常数流程

http://blog.sciencenet.cn/blog-567091-740196.html

• The phonon code ph.x calculates normal modes at a given q-vector, starting from data files produced by pw.x with a simple SCF calculation

• 4.2 Calculation of interatomic force constants in real space

First, dynamical matrices are calculated and saved for a suitable uniform grid of q-vectors

Second, code q2r.x reads the dynamical matrices produced in the preceding step and Fourier-transform them, writing a file of Interatomic Force Constants in real space, up to a distance that depends on the size of the grid of q-vectors

Program matdyn.x may be used to produce phonon modes and frequencies at any q using the Interatomic Force Constants file as input.

Since v.5.0, there are two ways of calculating electron-phonon coefficients, distinguished according to the value of variable electron_phonon.

The following holds for the case electron_phonon= 'interpolated' (see also Example 03).

It is convenient to use a coarse k-point grid to calculate phonons on a

suitable wavevector grid;

a dense k-point grid to calculate the sum at the Fermi energy.

The calculation proceeds in this way:

(1)  a scf calculation for the dense -point grid (or a scf calculation followed by a non-scf one on the dense -point grid); specify option la2f=.true. to pw.x in order to save a file with the eigenvalues on the dense k-point grid.

(2) a normal scf + phonon dispersion calculation on the coarse k-point grid, specifying option electron_phonon='interpolated', and the file name where the self-consistent first-order variation of the potential is to be stored: variablefildvscf).

(3) Finally, you can use matdyn.x and lambda.x (input documentation in the header of PHonon/PH/lambda.f90) to get the F( ) function, the electron-phonon coefficient , and an estimate of the critical temperature Tc

7)   计算超导温度流程

flfrc是用来设置保存所计算出来的力常数矩阵的文件名

fildyn是用来设置上一步已经计算得到的q点的动力学矩阵的文件名 Fig. S15. Phonon properties and Eliashberg spectral function for YH9 (a) and YH10 (b)

within the clathrate structures.

Phonon dispersion relations, projected phonon densities of states

(PDOS), and the Eliashberg spectral function (α2F(ω)/ω) and its integration λ(ω)【 λ取积分后的曲线走平后的稳定值？】 are plotted in the

left, middle, and right panels, respectively.

Circles in the phonon dispersion curves denote the

phonon linewidth with a radius proportional to the electron-phonon coupling strength.【声子带宽如何求，如何画图？】

Phonons

with large linewidths for YH10 at the Γ point correspond to the T2g mode at 1600 cm.1, the Eg

mode at 2000 cm.1, the T2g mode at 2200 cm.1, and the A1g mode at 2700 cm.1, respectively.【如何给出各个振动模式】

The

absence of imaginary frequencies in the phonon spectra indicates that all of the predicted clathrate

structures are dynamically stable.

example03:
This example shows how to calculate electron-phonon interaction
coefficients at X for fcc Al.

###Step1: SCF at dense k-mesh, good enough for electronic DOS,  specifying  la2F = .true.,  : pw.x

K_POINTS {automatic}
16 16 16  0 0 0

la2F = .true.,

###Step2: SCF at k-mesh good enough for phonons:  pw.x

K_POINTS {automatic}
8 8 8  0 0 0

####Step3:Specifying option electron_phonon='interpolated'. and the file name : fildvscf    : ph.x

cat > al.elph.in << EOF
Electron-phonon coefficients for Al
&inputph
tr2_ph=1.0d-10,
prefix='al',
fildvscf='aldv',
amass(1)=26.98,
outdir='$TMP_DIR/', fildyn='al.dyn', electron_phonon='interpolated', trans=.true., ldisp=.true. nq1=4, nq2=4, nq3=4 / EOF ##### Step4: Post processing: q2r and matdyn 后处理程序这么长，【每一步代表什么含义哈？】 # a) q2r and matdyn # cat > q2r.in << EOF &input zasr='simple', fildyn='al.dyn', flfrc='Al444.fc', la2F=.true. / EOF q2r.x < q2r.in > q2r.out cat > matdyn.in.freq << EOF &input asr='simple', amass(1)=26.98, flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.false. / 19 0.000 0.0 0.0 0.0 0.125 0.0 0.0 0.0 0.250 0.0 0.0 0.0 0.375 0.0 0.0 0.0 0.500 0.0 0.0 0.0 0.750 0.0 0.0 0.0 1.000 0.0 0.0 0.0 0.825 0.125 0.125 0.0 0.750 0.250 0.250 0.0 0.625 0.375 0.375 0.0 0.500 0.500 0.500 0.0 0.325 0.325 0.325 0.0 0.250 0.250 0.250 0.0 0.125 0.125 0.125 0.0 0.000 0.000 0.000 0.0 0.125 0.125 0.000 0.0 0.250 0.250 0.000 0.0 0.325 0.325 0.000 0.0 0.500 0.500 0.000 0.0 EOF matdyn.x < matdyn.in.freq > matdyn.out.freq cat > matdyn.in.dos << EOF &input asr='simple', amass(1)=26.98, flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.true. fldos='phonon.dos', nk1=10, nk2=10, nk3=10, ndos=50 / EOF matdyn.x < matdyn.in.dos > matdyn.out.dos #### b) Get Tc 【如何得到带权重的k点】 cat > lambda.in << EOF 10 0.12 1 ! emax (something more than highest phonon mode in THz), degauss, smearing method 8 ! Number of q-points for which EPC is calculated, 0.0000000 0.0000000 0.0000000 1.00 ! the first q-point, use kpoints.x program to calculate -0.2500000 -0.2500000 0.2500000 8.00 ! q-points and their weight -0.5000000 -0.5000000 0.5000000 4.00 ! 0.0000000 0.0000000 0.5000000 6.00 ! 4th q-point, qx,qy,qz -0.2500000 -0.2500000 0.7500000 24.00 ! -0.5000000 -0.5000000 1.0000000 12.00 ! 0.0000000 0.0000000 1.0000000 3.00 ! -0.5000000 0.0000000 1.0000000 6.00 ! the last q-point elph. 0.000000. 0.000000. 0.000000 ! elph output file names, elph.-0.250000. 0.250000.-0.250000 ! in the same order as the q-points before elph. 0.500000.-0.500000. 0.500000 elph. 0.000000. 0.500000. 0.000000 elph. 0.750000.-0.250000. 0.750000 elph. 0.500000. 0.000000. 0.500000 elph. 0.000000.-1.000000. 0.000000 elph.-0.500000.-1.000000. 0.000000 0.10 ! \mu the Coloumb coefficient in the modified ! Allen-Dynes formula for T_c (via \omega_log) EOF lambda.x < lambda.in > lambda.out cat lambda.out 中有Tc数据 8) k点或q点权重如何获得 a. 在布里渊区里取q 点，通过做scf 从scf.out 文件中获取 （已有权重）；选取的q 点最好不要少于50 个。 如： number of k points= 35 Methfessel-Paxton smearing, width (Ry)= 0.0300 cart. coord. in units 2pi/alat k( 1) = ( 0.0000000 0.0000000 0.0000000), wk = 0.0039062 k( 2) = ( 0.0000000 0.0000000 0.0397404), wk = 0.0234375 k( 3) = ( 0.0000000 0.0000000 0.0794808), wk = 0.0234375 k( 4) = ( 0.0000000 0.0000000 0.1192213), wk = 0.0234375 ........ http://blog.sciencenet.cn/blog-567091-725813.html : 自洽计算、kpoints测试和ecut测试 http://blog.sciencenet.cn/blog-567091-727329.htmlgamma点的测试很重要 b. 在步骤1 里面的q 点里选取一个q（不是Gamma）点做Kmesh2 的测定，通 过这一测定我们可以确定nscf2 里Kmesh2 以及nsig 的值。【确定K点及nsig值，测试脚本是...？，nsig代表什么？表中数值是能量吗？】 c. 选取了q 点，得到了kmesh2, 我们来计算（run_elph）每一个q 点的电声相互 作用常数λq（计算量很大，请不要着急） d. 通过这些q 点的电声相互作用常数λq 来计算总的λ（run_lambda 输入文件）， 可得到alpha2F.dat 为谱函数，lambda.dat 中有lambda 和声子频率对数平均值 第一行：3.6 0.1 1 emax, degaussq(THz), ngaussq ! emax (THz) : alpha2F is plotted from 0 to "emax" in "nex" steps ! degaussq (THz) : gaussian smearing for sum over q ! NB: not the same used in phonon! ! ngaussq : 0 for simple gaussian, 1 for Methfessel-Paxton etc.47 第二行：47 （下面的47 个q 点及权重） 0.0000000E+00 0.0000000E+00 0.0000000E+00 2.0000001E-03 ………… elph_1 (以下文件里是对应上面47 个q 点的声子和lambda(q)) elph_2 ………… elph_47 其中elph_1….elph_47 文件可由（la_*.dyn 和la_*.elph.out）通过lamd_in_get.sh 得到 给电声耦合作用文件编号——lamd_in_get.sh： #!/bin/sh for a in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 do awk ' BEGIN{r=1e20} /Diagonalizing the dynamical matrix/{r=NR+2} NR==r{print$4, $5,$6,10,3}

/omega$$.*THz.*cm-1/ { printf " %s ", 4 } /omega\([^ ].*THz.*cm-1/ { printf " %s ", 3 } END{printf "\n"} ' la_a.dyn > test1 awk '/Gaussian Broadening:/,/lambda\(3$$=/' la_$a.elph.out >test2 cat test1 test2 >elph_$a

done

5. 通过McMillan 方程来计算超导转变温度Tc。（用cal_tc.f 小程序）  新版PWSCF不需要改步骤

（1）输入w_log

（2）输入lambda

(3) 输入mustar

9）声子谱如何计算与绘图

The superconducting critical temperature, Tc, was calculated using the Quantum Espresso (QE)37 program package to obtain the dynamical matrix and the electron-phonon coupling (EPC) parameters.

In the QE calculations, the H pseudopotential, obtained from the QE pseudopotential library, was generated by the method of Trouiller-Martins38 with 1s1 valence configurations

and for Sc we employed an ultrasoft pseudopotential from the GBRV library, 39, 40 which consisted of the 3s23p64s23d1 valence configurations. Both pseudopotentials include the PBE generalized gradient approximation.

Plane-wave basis set cutoff energies were set to 80 Ry for all systems.

The Brillouin-zone sampling scheme of Methfessel-Paxton41 using a smearing of 0.03 Ry and dense k-point grids was used for all calculations carried out on the Sc-H phases.

Density functional perturbation theory42 as implemented in QE was employed for the phonon calculations.

The EPC parameter, λ, was calculated using a set of Gaussian broadenings in steps of 0.005 Ry from 0-0.300 Ry【为什么要采用一系列的高斯展宽来做】.

The broadenings for which λ was converged to within 0.05 were between 0.015 and 0.055 Ry for all structures【展宽or  λ收敛！】. Tc has been estimated using the Allen-Dynes modified McMillan equation43 ,

where wlog is the logarithmic average frequency 【对数平均频率】and

μ*, the Coulomb pseudopotential, was taken as 0.1.

In addition, Tc was calculated by numerically solving the Eliashberg equations based on the spectral function【谱函数？】, a2F(w), obtained with QE and using μ*=0.1.44

 Space group Pressure(GPa) k-point q-mesh McMillan Eliashberg Allen-Dynes P6/mmm-ScH2 300 16*16*16 4*4*4 4/- - - P63/mmc-ScH3 400 12*12*6 4*4*2 <1/- - -

 I4/mmm-ScH4 120 8*8*6 2*2*2 92 163 105 I4/mmm-ScH4 250 12*12*8 4*4*2 68 78 72

 Im-3m-ScH6 350 12*12*12 4*4*4 135 169 149

In this work, convergence tests of the kinetic energy cutoff as 60

Ry and 16  16  16 Monkhorst–Pack k-point sampling mesh29

was applied for the Brillouin zone integration, using a smearing

parameter of 50 mRy.

Ultraso? plane-wave pseudopotentials30

for Sc and H were applied to calculate the relative properties of

ScH2 and ScH3, and 3p63d14s2 for Sc as well as 1s1 for H were

used as electronic con?gurations to obtain a self-consistent

atomic charge density.

Phonon dispersion and electron–

phonon calculations were obtained by applying linear response

theory.

The dynamical matrices were calculated firstly at 10

wave (q) vectors using a 4  4  4 q grid【动力学矩阵，10个波矢量，4 4 4 网格】 in the irreducible

wedge of the Brillouin zone, and then the dynamical matrices

were Fourier-transformed to real space, accompanied by the

construction of force constant matrices, which were used to

obtain phonon frequencies at arbitrary points.

The Monkhorst–

Pack (MP) mesh was set as 12  12  12 and the Fermi–Dirac

smearing width was 0.02 Ry.

The zone-center phonon mode

Eigen displacements, which represent the infrared and Raman

activity of the cubic phase scandium hydrides【中心声子模式位移】 (ScH2 and ScH3)

were obtained using an isotropy code developed by Stokes and

Hatch.31

The energy band structures of fcc-ScH2 and ScH3 were

calculated in the148

irreducible fcc Brillouin zone (BZ) using the

generalized gradient approximation of the Perdew, Burke and

Ernzerhof (GGA-PBE) method.32,33 All of our work was carried

out using the Quantum-Espresso (QE) package.34

Phys. Rev. B. 2017, 96(9), 094513.

The electron-phonon coupling parameter (EPC) of the stable

compounds are calculated within the framework of linear

response theory using the QUANTUM-ESPRESSO code .

We employ the norm-conserving pseudopotentials

with 1s1 and

3d14s2 as valence electrons for H and Sc, respectively.

The

convergence tests concluded that suitable calculation parameters

are 80 Ry for kinetic energy cutoff, 16 × 16 × 16 k-point

sampling mesh, and 4 × 4 × 4 q mesh in the Brillouin zone.

Proc. Natl. Acad. Sci. U. S. A. 2017, 114(27), 6990-6995.

The superconducting Tc of these hydrides is expected to increase with hydrogen content of the material.1 Previous theoretical studies suggested sodalite-structured hydrides (e.g. CaH6 and YH6) should exhibit high Tc superconductivity. Electron-coupling calculations for LaH4 based BCS theory give a relatively small value (λ = 0.43). The estimated Tc is 10 K at 300 GPa with the typical choice of the Coulomb potential of μ*=0.1.

Using the McMillan equation, we calculated39 Tc using the spectral function (α2F(w)), again with μ*=0.1 . Here, λ is the first reciprocal moment of the spectral function, α2F(w), and ω(q) is the weight of a q point in the first BZ, and the EPC spectral function α2F(w) is expressed in terms of the phonon line width γqj arising from electron-phonon coupling. In this equation, Nf is the electronic density of electron states at the Fermi level. The linewidth γqj of a phonon mode j at wave vector q, arising from EPC is given by where the sum is over the first BZ, with ΩBZ as the volume of the BZ, and ξkn are the energies of bands measured with respect to the Fermi ξF level at point k. Here, gjkn,k+qm is the electron-phonon matrix element for scattering from an electron in band n at wave vector k state to band m at wave vector k+q via a phonon with wave vector q.

The calculated electron-coupling parameter for LaH8 at 300 GPa is relatively large (λ = 1.12), and the calculated Tc is 131 K assuming μ*=0.1 (Eq. 1).

That result encouraged us to calculate the Tc of LaH10. In this structure, the predicted H-H distance is 1.15 Å at these pressures, indicating weak H…H interactions

and no clear distinction between stretching and bending vibrations. As a result, all H vibrations effectively participate in the electron-phonon coupling process, which appears to enhance high superconductivity.

The calculated electron-phonon coupling is quite large (λ = 2.2). Tc was estimated from the spectral function (α2F(w)) by numerically solving the Eliashberg equations15,17 with μ* = 0.1. Coulomb repulsion is taken into account in terms of μ* scaled to a cutoff frequency. At 250 GPa, the estimated Tc of LaH10 is 280 K with μ*=0.1. The calculated Tc is found to decrease with increasing pressure.

1) PRB1975年Allen-Dynes文章摘录

The Eliashberg equations.

The McMillan equation  Following McMillan, we have calculated the electron-phonon coupling with the formula The  quantity which is often referred to as the Hopfield parameter, was calculated using the RMTA. This quantity, as shown in Fig.3, increases with increasing pressure, in contrast with N(EF), indicating a rapid increase of the electron-ion matrix element <I2>.

The matrix element <I2> is calculated from the angular momentum components of the DOS, the scattering phase shifts, and the free-scatterer DOS, all evaluated at EF. The detailed procedure followed is given in Ref. 7.

For the determination of the average phonon frequency we use the Debye-like approximation

the Allen-Dynes modified McMillan equation

The McMillan equation (with a prefactor altered from ΘD/1.45 to ωlog/1.2) is found to be highly accurate for all known materials withλ< 1.5 but in error for large values ofλ.

Correction factors to McMillan's equation are found in terms ofλ,μ*, and one additional parameter, (<ω2>)1/2/ωlog.

The frequency ωlog is defined as exp<lnω> where the averages<lnω> and 2> are defined using (2/λω)α2F(ω) as a weight factor.

We find several aspects of his work which agreed mell with information available at the time on medium-coupling superconductors (0. 5 <λ< 1)

but need modification for strongly coupled materials (λ>1) which have been more recently studied.

We attempt to make the appropriate modifications.

the electron-phonon spectral function α2F(ω) defined as The parameter λ is a dimensionless measure of the strength of α2F(ω)【电声耦合常数强度的量度，是标量】 http://blog.sciencenet.cn/blog-567091-1088020.html

全部精选博文导读

GMT+8, 2019-6-17 20:50