QUANTUM ESPRESSO DFPT EPW Wannier

QUANTUM ESPRESSO：EPW

Posted by Xiehua on December 16, 2021

QUANTUM ESPRESSO：epw.x计算electron-phonon coupling

EPW is an open-source community code for ab initio calculations of electron-phonon interactions using Density-Functional Perturbation Theory and Maximally Localized Wannier Functions.

计算步骤如下：

  1. 第一步：进入phonon目录进行scf自洽计算
    pw.x <scf.in> scf.out
    NOTE: outdir需要设置为”./”,由于后面需要对ph.x计算的结果使用pp.py进行收集。对于二维材料或者低维材料，建议在scf计算中使用参数：assume_isolated='2D'或者其他参数。

  2. 使用ph.x进行DFPT计算，得到dvscf文件：
    Note: set:

    1
2
3

     outdir = './'
 fildyn = 'graphene.dyn'
 fildvscf = 'dvscf'

    • 第二步：ph.x进行DFPT计算（最费时间，需要注意设置参数fildyn和fildvscf)
      在phonon计算中可以使用 -nimage N 参数进行并行计算，可以在不同的节点上计算不同的q点，也可以不同的任务提交不同的q点，然后进行后续处理，这样比较省机时，但是处理起来比较麻烦。On parallel machines the q point and the irreps calculations can be split automatically using the -nimage flag. See the phonon user guide for further information.

      1

      mpirun -np $NP -machinefile ${CURDIR}/nodelist ph.x -ni 6 -npool 28 <ph.in> ph.out`

  3. 第三步：使用pp.py收集ph.x计算得到的fildvscf相关文件到save文件夹。对于使用paw赝势计算的结果，需要采用修改版的pp-paw.py进行处理，或者使用新版本的pp.py。ph.x计算完成后使用q2r.x将力常数矩阵转换到实空间，并使用matdyn.x计算得到声子谱，检查声子谱是否正确，这是确保EPW计算没有问题的关键保证。

  4. 第四步进入epw目录，先进行scf计算（或者将phonon目录中的内容拷贝过来），再进行nscf计算，需要采用homogeneous kpoint grid进行nscf计算

    kmesh.pl 12 12 1 >>${prefix}.nscf.in

    需要修改calculation=nscf、nbnd=
    Note：nbnd的设置需要使得每个k点都包含所有的投影波函数的信息，通过fatband的能带图去看nbnd是否足够。

  5. 第五步设置epw.in文件，进行epw计算.先不用计算任何性质，主要是为了确定得到wannier函数求取成功，并得到wannier表象下的电神耦合矩阵，用于后续计算。注意fsthick的设置，会影响打印出来的电声耦合矩阵元包含的能带数和q点数。 note: dis_win_max的值要根据fatband的结果进行设置，不能随意设置的无限高或者不设置，这样会导致得到的wannier函数很难局域化。四个能量窗口都需要根据fatband的结果仔细的设置。

epw.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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

epw calculation of graphene
&inputepw
  prefix = 'graphene'
  outdir = './'
  amass(1)= 12.0107
  dvscf_dir = '../phonon/save'

  iverbosity = 0

  elph        = .true.
  epbwrite    = .false.
  epbread     = .false.
  epwwrite    = .true.
  epwread     = .false.
  etf_mem     = 1
  prtgkk   = .false.

!  eig_read    = .true.

!  asr_typ     = 'one-dim'
!  use_ws      = .true.

  wannierize = .true.
  nbndsub     =  6
!  bands_skipped = 'exclude_bands = 1-2'
  num_iter = 30000
  iprint   = 2
  dis_win_max = 25.0
  dis_win_min = -25.0
  dis_froz_min = -4.5
  dis_froz_max = 1.4
  proj(1) = 'C:sp2'
  wannier_plot= .true.
  wannier_plot_supercell = 5 5 1
!  wdata(1)= 'bands_plot = .true.'
!  wdata(2)= 'begin kpoint_path'
!  wdata(3)= 'G 0.00 0.00 0.00  M 0.50 0.00 0.00'
!  wdate(4)= 'M 0.50 0.00 0.00  K 0.3333 0.3333 0.00'
!  wdate(5)= 'K 0.3333 0.3333 0.00  G 0.00 0.00 0.00'
!  wdata(6)= 'end kpoint_path'
!  wdata(7)= 'bands_plot_format = gnuplot'
  wdata(8)= 'conv_tol      = 1.0e-10 '
  wdata(9)= 'conv_window   = 10      '
  wdata(10)= 'dis_conv_tol  = 1.0e-10 '
  wdata(11)= 'dis_conv_window = 5     '
  wdata(12)= 'dis_num_iter= 30000      '
  wdata(13)= 'dis_mix_ratio= 0.5      '
  wdata(14)= 'guiding_centres = .true.'
  wdata(15)= 'translate_home_cell  : true'
  wdata(16)= 'translation_centre_frac :   0.0 0.0 0.0  '

  elecselfen  = .false.
  phonselfen  = .false.
  a2f         = .false.

  fsthick     = 5.0 ! eV
  temps       = 300 ! K
  degaussw    = 0.005 ! eV

!  band_plot   = .true.
!  filkf       = './LGX.txt'
!  filqf       = './LGX.txt'

  nkf1 = 12
  nkf2 = 12
  nkf3 = 1
  nqf1 = 12
  nqf2 = 12
  nqf3 = 1

  nk1 = 12
  nk2 = 12
  nk3 = 1
  nq1 = 12
  nq2 = 12
  nq3 = 1
/

作业提交bsub脚本：qe-epw.bsub

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

#!/bin/bash
#BSUB -J epw
#BSUB -q privateq-zw
##BSUB -q publicq
#BSUB -n 28
#BSUB -R "span[ptile=28]"
#BSUB -o %J.out
#BSUB -e %J.err

CURDIR=$PWD
rm -f nodelist >& /dev/null
for host in `echo $LSB_HOSTS`
do
echo $host >> nodelist
done
NP=`cat nodelist | wc -l`

prefix='graphene'
#source ~/xiehua/.bashrc
#export OMP_NUM_THREADS=1
#export MKL_NUM_THREADS=1
export MODULEPATH=/share/home/zw/xiehua/modulefiles:$MODULEPATH

module load Quantum_Espresso/7.0

#mpirun -np $NP pw.x -nk 7 <vc-relax.in>vc-relax.out
#mpirun -np $NP pw.x -nk 7 <scf.in>scf.out
#mpirun -np $NP pw.x -nk 7 <nscf.in> nscf.out
#wannier90.x -pp ${prefix}
#pw2wannier90.x < pw2wan.in > pw2wan.out
#mpirun -np $NP wannier90.x ${prefix}
mpirun -np $NP epw.x -npool $NP <epw.in > epw.out

  • EPW声子谱和QE不一致 (参考简书)
    这个问题一般是由于声子求和规则导致的，EPW中提供了读入实空间力常数来计算声子频率的方法，并且也提供了相应的声子求和规则（与matdyn.f90里面的相同）。只需要改lifc = .true.，然后再设置声子求和规则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,在q2r.x的计算中，需要设置参数flfrc='prefix.fc',这样在使用pp.py进行处理时，会自动将prefix.fc保存进save目录，并改名为ifc.q2r 强烈建议：在EPW计算中，使用ph.x的结果后，使用q2r.x计算的结果进行EPW计算，这样可以保证满足声学支求和规则和LO-TO分裂。

    1
2
3

      &input  
  fildyn='carbyne.dyn', zasr='simple', flfrc='carbyne.fc',loto_2d=.true.  
  /

6.采用5步中计算成功后的结果，使用第五步epwwrite=.true.设置输出的wannier表象下的电声耦合文件，计算不同nqf和nkf插值密度的EPW计算。采用脚本runepw.sh计算其他插值q点和k点情况下的电神耦合矩阵元，并输出gmnvkq矩阵元。

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

#!/bin/bash
for i in $(seq 6 6 36)
  do
  cp epw.in epw${i}.in
  sed -i "s:epwwrite:epwwrite=.false. ! :g" epw${i}.in
  sed -i "s:epwread:epwread=.true. !:g" epw${i}.in
  sed -i "s:epbread:epbread=.false. !:g" epw${i}.in
  sed -i "s:epbwrite:epbwrite=.false. !:g" epw${i}.in
  sed -i "s:prtgkk:prtgkk=.true. !:g" epw${i}.in
  sed -i "s:wannierize:wannierize=.false. !:g" epw${i}.in
  sed -i "s:wannier_plot:wannier_plot=.false. !:g" epw${i}.in
  sed -i "s:elecselfen:elecselfen=.false. !:g" epw${i}.in
  sed -i "s:fsthick:! fsthick:g" epw${i}.in

  sed -i "s:nkf1:nkf1=$i !:g" epw${i}.in
  sed -i "s:nkf2:nkf2=$i !:g" epw${i}.in  
# sed -i "s:nkf3:nkf3=$i !:g" epw${i}.in   
  sed -i "s:nqf1:nqf1=$i !:g" epw${i}.in
  sed -i "s:nqf2:nqf2=$i !:g" epw${i}.in  
# sed -i "s:nqf3:nqf3=$i !:g" epw${i}.in 

  cp qe-epw.bsub qe-epw${i}.bsub
  sed -i "2s:epw:epw${i}:g" qe-epw${i}.bsub
  sed -i "s:epw.in:epw${i}.in:g" qe-epw${i}.bsub
  sed -i "s:epw.out:epw${i}.out:g" qe-epw${i}.bsub	

  bsub < qe-epw${i}.bsub
  sleep 10
  done
  1. 采用较密的k点和q点计算电子自能，用于计算电子弛豫时间。使用脚本runepwmkq.sh
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

#!/bin/bash
for i in $(60)
  do
  cp epw.in epw${i}.in
  sed -i "s:epwwrite:epwwrite=.false. ! :g" epw${i}.in
  sed -i "s:epwread:epwread=.true. !:g" epw${i}.in
  sed -i "s:epbread:epbread=.false. !:g" epw${i}.in
  sed -i "s:epbwrite:epbwrite=.false. !:g" epw${i}.in
  sed -i "s:prtgkk:prtgkk=.false. !:g" epw${i}.in
  sed -i "s:wannierize:wannierize=.false. !:g" epw${i}.in
  sed -i "s:wannier_plot:wannier_plot=.false. !:g" epw${i}.in
  sed -i "s:elecselfen:elecselfen=.true. !:g" epw${i}.in
  
  sed -i "s:nkf1:nkf1=$i !:g" epw${i}.in
  sed -i "s:nkf2:nkf2=$i !:g" epw${i}.in  
# sed -i "s:nkf3:nkf3=$i !:g" epw${i}.in   
  sed -i "s:nqf1:nqf1=$i !:g" epw${i}.in
  sed -i "s:nqf2:nqf2=$i !:g" epw${i}.in  
# sed -i "s:nqf3:nqf3=$i !:g" epw${i}.in 

  cp qe-epw.bsub qe-epw${i}.bsub
  sed -i "2s:epw:epw${i}:g" qe-epw${i}.bsub
  sed -i "s:epw.in:epw${i}.in:g" qe-epw${i}.bsub
  sed -i "s:epw.out:epw${i}.out:g" qe-epw${i}.bsub	

  bsub < qe-epw${i}.bsub
  done

8.EPW运算过程中可能遇到的报错问题。
（1）在epw.x运算过程中可能遇到报错如下：

1
2
3
4
5
6
7
8
9

     Fermi energy coarse grid =  -1.703677 eV

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine efermig (1):
     internal error, cannot bracket Ef
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...

该问题可能导致的原因，参考EPW Forum 可能的解决方案包括：1.使用更多的Wannier函数进行求解局域的WFs，拟合的wannier能带加上exclude_bands需要包含费米面一下的所有能带。2.使用bands_skipped = 'exclude_bands = 1-16'排除部分能带。3.dis_win_min一般不需要设置。4.使用的wannier函数数目应该大于等于体系中的电子数目，否者在计算Ef时会有错误，这是由于EPW在计算占据态电子数目的bug导致的。问题出在计算插值后的费米能量时，nelec设置有误，需要将nelec修改为wannier函数中的占据态电子数，nelec_wf,或者从输入文件中读入费米能，最简单的方法就是在epw.in中根据nscf计算的结果指定费米能。 在epw.in中设置参数efermi_read=.true.和fermi_energy=，但是这样如果使用剪刀差时，epw计算得到的价带顶和导带底的带指标是错误的。最好的方法是，在Wannier拟合的过程中，必须包含费米面一下的所有能带，对于不能包含的能带，必须使用bands_skipped='exclude_bands='进行排除。

例如，对于graphene的epw的计算，需要设置proj(1) = 'C:pz,px,py',虽然Wannier拟合只需要采用proj(1) = 'C:pz'即可完成拟合。

在epw计算中，其计算得到的体系的电子数可能是一个错误的，导致计算出来的Ef也是一个错误的，这是epw的一个bug。 pew forum Currently, since your nbndsub (# of Wannier functions) is smaller than the number of electrons in your system, you encounter the issue with calculation of Fermi energy on fine grids.

nbndsub is the number of Wannier functions and you need to determine it based on the band manifold of interest. To determine the Fermi energy, nbndsub should be equal to or larger than the number of valence electrons.

EPW的bug： EPW计算的Ef有问题：只有在wannier拟合的第一条band正好是Bloch空间的第一条band才是正确结果。可以选择将不需要拟合的能带排除出去，或者将低能带的轨道全部拟合。该报错只出现在EPW中，在Wannier90中不存在这个问题。

EPW中计算，必须保证Wannier拟合的能带加上exclude_bands的能带包含费米面一下的所有能带，否则在计算插值后的费米能时会有bug。

CATALOG
    FEATURED TAGS
    LVCSH QUANTUM ESPRESSO Wannier Fortran QE 激活秘钥 DFPT 安装 EPW Endnote X20 Visual Studio WSL phonon
    FRIENDS