电声耦合计算

Quantum Espresso PHonon模块

GRID example

网格的大小由下面几个参数决定
nqs q点的数量
x_q q点的坐标
nfs 频率（imaginary frequencies）的个数
fiu 哪一个频率
在recover计算中，程序首先从phq_readin读入文件（如果这是起始的计算，那phq_readin会设置这些参量）frequency和tensors从输入文件中读。
comp_iq 决定这个q点是否计算
comp_irr_iq决定这个表示是否计算
comp_iu这个频率是否计算

EPW

degaussw
电子在晶格中运动的时候会使晶格发生微小的畸变（由于库伦相互作用），晶格畸变反过来又作用到电子上，是的电子动力学发生变化，导致电子的quasi-particle state有效质量增加，lifetime降低，在场论里面，人们用电子自能来描述由于电声耦合导致的电子动力学的变化。自能的实部描述的是电子能量的变化，虚部描述的是电子lifetime的变化

电子自能算符的对角项作用在unperturbed 电子态上得到

计算电子自能的时候公式里边

Adapted from slides of S. Ponre

Beta相硼烯测试结果.png

Chi相硼烯测试结果.png

eptemp is the temperature at which you want your results. The temperature eptemp will give the value of the Bose-Einstein n(T) factors.

degaussw is the broadening on the delta of electron properties.
degaussq is the broadening on the delta of phonon properties.
If those parameters are too large, you will wash out important features. If they are too small you will have numerical instabilities.
In theory, those two parameters can be 0 if you have dense enough k/q-points grids.

In my opinion, the correct way to converge them is to do the following:

1. set degaussw and degaussq to relatively large value and then increase the k/q sampling until your results converge
2. decrease degaussw and degaussq and then increase the k/q sampling until your results converge (you will notice that you need larger k/q grid >to achieve convergence).
3. repeat until the difference of the converge results (in k/q) converge also in degaussw and degaussq.

Note that if you study for example Eliashberg spectral functions, different value of degaussq will be outputed by default in a2F file.
I usually converge on degaussw only and check that degaussq is relatively stable.

delta_approx
这个标签决定在计算声子线宽(linewidth)的时候是否采用double delta approximation
根据EPW forum上roxana(EPW tutorial里面讲超导的主讲人)描述，EPW程序里面关于超导性质的估算和能隙的计算都是采用的double delta approximation的。
我在计算中测试了使用delta_approx为true和false的两种情况，在k点足够密集以后都得到了收敛的结果，但是两者明显不同。(y?)

1.jpg

输出文件及其包含的内容

prefix.epb 包含最初的k/q点网格上的Hamiltonian、动力学矩阵元、电声耦合矩阵元。
prefix.epmatwp1，crystal.fmt，dmedata.fmt，epwdata.fmt 转化到Wannier表象中的Hamiltonian、动力学矩阵元、电声耦合矩阵元。
ep_coupling
elph
这两个参数是用来计算*.ephmat, *.freq, *.egnv, and *.ikmap这些文件的，如果前面已经得到了这些文件，那么在接下来计算超导的时候就可以把这两个参数关掉，但是用到的cpu核数必须和之前一致，因为ephmat文件个数和核数一样。具体参见EPW-forum
epbread这个是读取bloch表象的电声耦合矩阵元的标签，实际上有了 *.epmatwp1文件之后，就不需要读取epb了，而是直接读取 *.epmatwp1和epwdata.fmt文件，并且读取 *.epmatwp1文件没有核数的限制，它是wannier表象的电声耦合矩阵元，设置kmaps = .true.就可以了，但是这里其实并没有读取"prefix.kmap" and "prefix.kgmap"这两个文件。只是读取了 *.epmatwp1和epwdata.fmt文件。
ephwrite
这个参数是用来控制是否输出 *.ephmat文件的，这个文件里包含了在Fermi window里fine k、q mesh上前面用 elph,ep_coupling计算出来的电声耦合矩阵元， *.ephmat文件个数和使用的核的个数相同，这个文件和 *.freq、 *.egnv（分别包含Fermi window里面的声子和电子本征值） *.ikmap（包含Fermi window里面的不可约k点的坐标）加在一起这四个文件包含了求解anisotropic Eliashberg方程的所有信息，求解其他温度的AE方程的时候也会用到这几个文件，但是如果你改了fsthick或者k、q点网格或者是使用cpu核的个数的时候这些将无法reuse。

读入Eliashberg谱函数求各向同性Eliashberg方程

fila2f = 'prefix.a2f'
EPW中提供了一种直接通过Eliashberg谱函数求解各向同性Eliashberg方程的方法，只需要提供的信息就行了，不过文件的格式以及单位必须和EPW自己产生的文件一致，第一列是声子频率，单位是meV，第二列是谱函数，应该是无量纲数。同时注意控制读写的输入参数应该与这个里面一致 Pade approximation

&inputepw
prefix = 'BaTiO3',
amass(1) = 137.34
amass(2) = 47.90
amass(3) = 16.00
outdir = './'

ep_coupling = .false.
elph = .false.
kmaps = .false.
epbwrite = .false.
epbread = .true.
epwwrite = .false.
epwread = .true.


fsthick = 0.5 ! eV
eptemp = 0.075 ! K
degaussw = 0.05 ! eV

a2f = .false.

dvscf_dir = '../phonons/save'

fila2f = 'BaTiO3.a2f'

ephwrite = .false.
eliashberg = .true.

laniso = .false.
liso = .true.
lreal = .false.
limag = .true.
lpade = .true.

conv_thr_iaxis = 1.0d-3
wscut = 0.525 ! eV

nstemp = 2
tempsmin = 0.005
tempsmax = 0.030

nsiter = 500

muc = 0.1

nkf1 = 20
nkf2 = 20
nkf3 = 20

nqf1 = 20
nqf2 = 20
nqf3 = 20

nk1 = 4
nk2 = 4
nk3 = 4
/
10 cartesian
0.0000000 0.0000000 0.0000000 0.0312500
0.0000000 0.0000000 0.2500000 0.1875000
0.0000000 0.0000000 -0.5000000 0.0937500
0.0000000 0.2500000 0.2500000 0.3750000
0.0000000 0.2500000 -0.5000000 0.3750000
0.0000000 -0.5000000 -0.5000000 0.0937500
0.2500000 0.2500000 0.2500000 0.2500000
0.2500000 0.2500000 -0.5000000 0.3750000
0.2500000 -0.5000000 -0.5000000 0.1875000
-0.5000000 -0.5000000 -0.5000000 0.0312500

EPW声子谱和QE不一致

这个问题一般是由于声子求和规则导致的，EPW中提供了读入实空间力常数来计算声子频率的方法，并且也提供了相应的声子求和规则（与matdyn.f90里面的相同）。只需要改lifc = .t.，然后再设置声子求和规则asr_typ = crystal（我一般都取crystal），同时需要注意的是要保证之前计算QE得到的文件通过pp.py收集起来那个必须有q2r.x产生的实空间力常数文件并且已经被命名为ifc.q2r，对于包含SOC的情况，这个文件必须叫ifc.q2r.xml并且是xml格式的文件。（这个一般不是太老的脚本pp.py都会自动帮你做这件事情。）参考phonon bandstructure from EPW and matdyn.x don't match

电声耦合常数偏低的问题

这段时间被这个问题所困扰，无法重复出文献中的数值，doping之后的单层电声耦合计算总是偏低，后来发现是smearing和层间距的问题，这里简单介绍一下电子结构计算中smearing的选取。
首先我们要明确为什么要有电子展宽，对于DFT里面很多参量（total energy、charge density）的计算，需要对占据态做求和，求和的过程中就会发现如果我们按照严格的基态的Fermi-Dirac分布来看，费米面以上的占据数严格为0的话，那么我们往往需要非常密集的K点sampling才能取得收敛的结果，因为费米面附近的精度将会大大影响计算结果，是否计入某个点可能会使结果变化很大。为了克服这一点，人们提出使用展宽的方式来使得我们即使在不那么密集的k点取样的情况下也能得到和严格情形下密集取样类似的结果。详见theos-ElectronicTemperature

tetrahedron & tetrahedron method with Blochl correction
这个方法适合计算体相材料的总能和态密度，这个方法没办法做分数占据。所以计算金属的原子受力和压力张量会有5%到10%的偏差。
mp
可以说是对高斯展宽的一般化，0阶mp分布就对应于高斯展宽。这个方法适合声子的计算，对这种方法的简单介绍可以参考这篇文章Methfessel-Paxton，大概思想就是用高阶厄密多项式来展开费米面附近的展宽，这个方法也可以对总能有很好的估计，但是展宽值的选取需要格外小心，展宽太大算出来总能可能不准确，小的展宽需要比较密集的k点取样。一个参考标准就是自由能与总能之差小于1meV/atom(针对VASP中的计算）。对于比较大的超胞MP方法也是很好的选择。不适用于半导体和绝缘体。详见VASP-ISMEAR，需要注意的是，这里虽然取的展宽，但是计算的总能是基态的总能也就是对应的0温的总能。但是mp展宽有时候会出现负占据和大于1占据的问题。

marzari-vanderbilt
Marzari为了解决上面的负占据和大于1占据的问题构造出来的方法，也叫做cold-smearing。我在QE的计算能带的example里面看到经常使用mv展宽，但是最近就是在这个上面不小心导致计算doping 的声子的软化和电声耦合远小于文献中的数值，所以计算声子的时候还是尽量使用mp展宽。
Fermi-Dirac
这个按照道理来说是最接近有限温情形下的电子分布的，但是使用这个也有一些问题，比如说想要得到比较收敛的结果需要比较大的展宽（0.1~0.5eV），这时候Fermi-Dirac分布的尾巴就会比较长，就需要算入很多的态，增加计算量。所以也不是说这个就比别的好，有时候可能还不如用一个比较假的smearing比如高斯。