使用QEtoolkit便利地创建QE的各种输入文件

丁越

使用QEtoolkit便利地创建QE的各种输入文件

------------------------------------------------------
2022.9.11 注：修正了一些问题
2022.8.10 注： 修正了DFT+U部分的问题
-----------------------------------------------------

一、前言

  Multiwfn程序中强大的输入文件生成功能给用户带来了极大的便利，特别是CP2K输入文件生成的功能，成为了该用户群体的一大福音。相对于CP2K来说，虽然QE的输入文件不难写，里面包含的内容也比较简明，但是如果能将输入文件生成进一步“傻瓜化”，避免在pwgui中点来点去落下某些关键参数的问题，这也将是一件非常有意义的工作。受此启发下，于是决定也搞一个QE的输入文件生成工具，继QEtoolkit-1.0版本后继续扩充成为了2.0版本，里面包含了pw.x，neb.x，dos.x，projwfc.x，hp.x，ph.x，pp.x，bands.x模块常用功能的输入文件生成。该工具的使用习惯和Multiwfn十分相似，所以想要生成什么输入文件按照提示输入相应选项、文件即可。第三节将会使用几个简单案例进行展示。

二、程序使用说明

  由于QE的模块比较杂，为了能够明白什么时候该输入什么文件，有必要先简单说明一下各个模块的用途：

  * pw.x: 这个是QE的核心程序，里面包含单点计算、非自洽计算、结构优化、能带计算、MD功能。

  * neb.x: 这个是利用CI-NEB计算过渡态的程序，neb.x本质上还是一个外壳，核心的自洽计算还是交给pw.x执行。

  * dos.x, projwfc.x: 这两者是用于DOS、PDOS计算的程序。

  * hp.x: 通过线性响应计算hubbard U的程序，当然还可以计算DFT+U+V。

  * ph.x: hp和ph乍一看挺混淆的，但是要知道这两个名字咋来的就记住了。hp刚才上面已经说了，是hubbard U。而这里ph是phonon，利用线性响应进行声子计算的程序。

  * pp.x: 该程序是用于后处理，顾名思义，post processing的程序。用于各种电子结构性质的计算，比如电子密度、静电势、RDG等等。

  * bands.x: 这个是pw.x 能带计算后进行后处理的程序，能生成gnuplot绘图格式的能带数据。

    在使用该工具包之前，用户首先需要更改QEtoolkit代码首行中user defined parameters, 包括默认赝势库的路径以及QE的版本。QEtoolkit预定义了四套赝势库，分别是GBRV的超软赝势、PSlibrary中pbe和pbesol的PAW赝势，ONCV的模守恒(NC)赝势，用户也应该自行去把这四套赝势备齐，这就足以应付绝大多数的计算问题了。一般计算首选GBRV的超软赝势，精度和速度都很不错，但是存在的问题是赝势并没有像PSlibrary中那么齐全，尤其是一些镧系元素上。其次PSlibrary中的赝势也是较好的选择，pbesol对于固体晶格参数的计算精度比较好。内置的PSlibrary中的赝势全部采用了非半芯、无非线性核校正赝势，如果你想用那些的话可以自行把预定义的赝势换了。带个ONCV的模守恒赝势主要是为了应付拉曼光谱计算，QE比较坑爹的一点就是好多计算没给你说明不能用某些泛函或者参数，比如这里的拉曼光谱计算必须就得用lda+NC去做。PS：也可以修改源码采用pbe+NC去做（https://minerals.chazeon.com/post/notes-on-raman-with-qe/）。QE的版本参数需要定义的原因是自7.1版本起DFT+U关键词进行了完全的重构，所以不同版本这块关键词不能混用。

三、案例

2.1 BCC Fe晶胞优化

  直接输入./QEtoolkit.sh Fe.cif，或者先运行QEtoolkit，选择功能0后再按照提示输入Fe.cif。接下来输入下面选项来生成晶胞优化的输入文件：

1. 0 #选择pw.x模块
   
2. 1 #选择计算体系的类型
   
3. 3 #这里Fe是导体
   
4. 2 #选择计算任务
   
5. 4 #晶胞优化
   
6. 6 #进入定义磁矩选项，弹出以下界面
   
7. Define Magnetization:
   
8. type1: Fe      magnetization:0         Natms:2
   
9. 
   
10. * Magnetization is defined as (nalpha-nbeta)/(nalpha+nbeta), ranges from -1 to 1.
    
11. * If you want to define antiferromagnetism state, two groups of same type atoms should be setted with opposite magnetization of -n and n, respectively.
    
12. * Input "q" to exit.
    
13. Which type of atoms do you want to define magnetization:
    
14. 1
    
15. Input magnetization, e.g. 0.5
    
16. 0.5
    
17. 
    
18. type1: Fe      magnetization:0.5       Natms:2
    
19. Which type of atoms do you want to define magnetization:
    
20. q
    
21. 这里只有Fe一种元素，定义完磁矩值后输入"q"退出定义磁矩界面
    
22. 8 #设置k点，这里定义3*3*3

复制代码

最后当前目录下就出现了Fe.vcrelax.in的输入文件，用户需要进一步修改截断能参数。当然这里直接使用了GBRV的赝势，你可以进入10选项来修改想要使用的赝势库类型。

2.2 Cu(111) 的功函数

  功函数公式为局部真空区域Hartree势减去Fermi能级。所以第一步先进行自洽计算得到Cu(111)面的Fermi能级，第二步通过pp.x产生静电势的Cube文件来获取局部真空区域Hartree势。
  第一步：生成Cu(111)的 单点计算文件

1. ./QEtoolkit Cu111.cif
   
2. 0 #pw.x 模块
   
3. 1 #Cu为导体
   
4. 3
   
5. 8 #设置k点，真空层方向保留一个k点即可
   
6. 3,3,1
   
7. 0 #生成输入文件

复制代码

运行结束后，输入grep Fermi Cu111.scf.out 即可获取Fermi能级为2.0330 eV。
  第二步：生成计算静电势的输入文件

1. ./QEtoolkit
   
2. 6 #进入pp.x模块
   
3. 5 #静电势的cube文件

复制代码

此时输入文件就被输出到了当前目录下，即esp.pp.in。用户需要手动修改里面的prefix为Cu111, 这必须与pw.x中的prefix保持一致，否则程序无法运行。此外，pp.x的运行命令也被输出到了屏幕上。运行结束后，我们将得到的esp.cube放到Multiwfn进一步处理。

1. Multiwfn.exe esp.cube
   
2. 13 #进入格点数据处理功能
   
3. 18 #计算沿着Z方向的面平均
   
4. Z
   
5. a
   
6. 3 #绘图
   
7. 6 #保存绘制的图
   
8. 9 #输出面平均的原始数据

复制代码

  这里不要直接目测读数据，我们看到图中12埃附近Hartree势就比较平直了，所以打开刚才输出的原始数据找到12埃附近的Z值，大约为0.51a.u. ,注意QE中的基本单位为里德堡(Ry)制，需要将这里得到的Z值转化为eV再去算功函数，所以功函数为0.51*13.6-2.0330=4.903 eV。

2.3 Si 的单点及能带计算

  第一步： Si晶胞的自洽计算

1. ./QEtoolkit Si.cif
   
2. 0
   
3. 1
   
4. 2 #Si为半导体
   
5. 8
   
6. 4,4,4
   
7. 0

复制代码

第二步：能带计算。首先生成能带计算的输入文件。

1. ./QEtoolkit Si.cif
   
2. 0
   
3. 1
   
4. 2 #设置计算任务类型为能带计算
   
5. 6
   
6. 0
   

复制代码

紧接着我们将上步产生的Si.scf.in拖到Seek-Path中获取高对称点路径，然后修改Si.bans.in中最后一行为获取到的路径。能带计算完成后，Si.bands.out中就有各个高对称点的x坐标，并且产生的bands.dat.gnu是适合gnuplot绘图的数据格式。我们再次启动QEtoolkit，在其他功能中的19功能，选项4是一个绘制band的脚本。如果是需要绘制其他结构的能带那就恰当修改里面脚本参数即可。

2.4 FeO的Hubbard U值计算   
  FeO是反铁磁绝缘体，如果不加U的话结果计算下来是导体。当前版本的QEtoolkit中无法直接定义反铁磁态，下个小版本更新时将会加入这个功能，所以当下只能先定义铁磁态，然后手动修改成反铁磁态。

1. ./QEtoolkit FeO
   
2. 0
   
3. 6
   
4. Define Magnetization:
   
5. type1: Fe      magnetization:0         Natms:2
   
6. type2: O       magnetization:0         Natms:2
   
7. 
   
8. * Magnetization is defined as (nalpha-nbeta)/(nalpha+nbeta), ranges from -1 to 1.
   
9. * If you want to define antiferromagnetism state, two groups of same type atoms should be setted with opposite magnetization of -n and n, respectively.
   
10. * Input "q" to exit.
    
11. Which type of atoms do you want to define magnetization:
    
12. 1
    
13. Input magnetization, e.g. 0.5
    
14. 0.5
    
15. 
    
16. type1: Fe      magnetization:0.5       Natms:2
    
17. type2: O       magnetization:0         Natms:2
    
18. Which type of atoms do you want to define magnetization:
    
19. q
    
20. 7 #开启DFT+U
    
21. 2
    
22. 8
    
23. 3,3,3
    
24. 0

复制代码

如上面提示所说，定义反铁磁态就要把Fe划分为两种，Fe1和Fe2，所以元素类型数目ntyp也由2改为3。U值也要设置一个非零值以需后续U值计算。最终更改完的输入文件如下所示：

1. &CONTROL
   
2.    calculation     = 'scf'
   
3.    restart_mode    = 'from_scratch'
   
4.    outdir          = './tmp'
   
5.    pseudo_dir      = '/home/dy/qe_pbe_GBRV'
   
6.    prefix          = 'FeO'
   
7.    verbosity       = 'low'
   
8. /
   
9. &SYSTEM
   
10.    ibrav           = 0
    
11.    nat             = 4
    
12.    ntyp            = 3
    
13.    ecutwfc         = 20.0
    
14.    ecutrho         = 300.0
    
15.    occupations     = 'smearing'
    
16.    degauss         = 0.001
    
17.    smearing        = 'gaussian'
    
18.    nspin           = 2
    
19.    starting_magnetization(1) = 0.5
    
20.    starting_magnetization(2) = -0.5
    
21.    lda_plus_u      = .true.
    
22.    Hubbard_U(1)    = 1.D-8
    
23.    Hubbard_U(2)    = 1.D-8
    
24.    U_projection_type = 'ortho-atomic'
    
25. /
    
26. &ELECTRONS
    
27.    electron_maxstep = 128
    
28.    conv_thr        = 1.D-6
    
29.    mixing_mode     = 'plain'
    
30.    mixing_beta     = 0.7
    
31.    mixing_ndim     = 8
    
32.    diagonalization = 'david'
    
33. /
    
34. CELL_PARAMETERS angstrom
    
35.      3.102293760    0.000000000    0.000000000
    
36.      1.499036090    2.721409730    0.000000000
    
37.      1.504046230    0.833573050    5.033304410
    
38. ATOMIC_SPECIES
    
39.    Fe1   55.845  fe_pbe_v1.5.uspp.F.UPF
    
40.    Fe2  55.845  fe_pbe_v1.5.uspp.F.UPF
    
41.    O    15.999  o_pbe_v1.2.uspp.F.UPF
    
42. ATOMIC_POSITIONS angstrom
    
43.    Fe1      0.000000000    0.000000000    0.000000000
    
44.    Fe2      3.052688000    1.777491000    2.516652000
    
45.     O      4.676533000    2.712576000    3.766550000
    
46.     O      1.428843000    0.842407000    1.266754000
    
47. K_POINTS automatic
    
48. 3 3 3 0 0 0

复制代码

启动QEtoolkit，选择3功能，然后选择1仅仅对一种元素进行响应计算U值。修改生成的hp.in中的prefix，其要与pw.x中的prefix保持一致。hp.x怎么运行屏幕上输出的提示就写的很清楚了，并且也详细输出了反铁磁绝缘体计算U值步骤。所以下面按照提示首先将FeO.scf.in复制为FeO.scf2.in，修改成如下所示：

1. &CONTROL
   
2.    calculation     = 'scf'
   
3.    restart_mode    = 'from_scratch'
   
4.    outdir          = './tmp'
   
5.    pseudo_dir      = '/home/dy/qe_pbe_GBRV'
   
6.    prefix          = 'FeO'
   
7.    verbosity       = 'low'
   
8. /
   
9. &SYSTEM
   
10.    ibrav           = 0
    
11.    nat             = 4
    
12.    ntyp            = 3
    
13.    ecutwfc         = 20.0
    
14.    ecutrho         = 300.0
    
15.    occupations     = 'fixed'
    
16.    tot_magnetization = 0
    
17.    nbnd = 26
    
18.    nspin           = 2
    
19.    starting_magnetization(1) = 0.5
    
20.    starting_magnetization(2) = -0.5
    
21.    lda_plus_u      = .true.
    
22.    Hubbard_U(1)    = 1.D-8
    
23.    Hubbard_U(2)    = 1.D-8
    
24.    U_projection_type = 'ortho-atomic'
    
25. /
    
26. &ELECTRONS
    
27.    electron_maxstep = 128
    
28.    conv_thr        = 1.D-6
    
29.    mixing_mode     = 'plain'
    
30.    mixing_beta     = 0.7
    
31.    mixing_ndim     = 8
    
32.    diagonalization = 'david'
    
33.    startingpot = 'file'
    
34.    startingwfc = 'file'
    
35. /
    
36. CELL_PARAMETERS angstrom
    
37.      3.102293760    0.000000000    0.000000000
    
38.      1.499036090    2.721409730    0.000000000
    
39.      1.504046230    0.833573050    5.033304410
    
40. ATOMIC_SPECIES
    
41.    Fe1   55.845  fe_pbe_v1.5.uspp.F.UPF
    
42.    Fe2  55.845  fe_pbe_v1.5.uspp.F.UPF
    
43.    O    15.999  o_pbe_v1.2.uspp.F.UPF
    
44. ATOMIC_POSITIONS angstrom
    
45.    Fe1      0.000000000    0.000000000    0.000000000
    
46.    Fe2      3.052688000    1.777491000    2.516652000
    
47.     O      4.676533000    2.712576000    3.766550000
    
48.     O      1.428843000    0.842407000    1.266754000
    
49. K_POINTS automatic
    
50. 3 3 3 0 0 0
    

复制代码

计算完后接下来就是最后一步hp.x进行线性响应计算U值了。hp.x计算结束后，就输出了名为FeO.Hubbard_parameters.dat的文件，里面Fe1和Fe2的U值分别为5.58和5.97.

2.5 TiO2(110)-H2O的振动频率及其热力学性质计算
  由于QE的声子计算不支持DFT-D3和DFT-D3(BJ)，所以为了后续的声子计算，改用了在DFT-D2下进行结构优化，这一点也确实让人难受，这都什么年头了还让人用DFT-D2，不知道下个版本会不会修正这个问题。生成结构优化文件的过程这里就不再赘述了，下面我们直接生成gamma点声子频率计算的输入文件，。

1. ./QEtoolkit.sh TiO2_H2O.relax.in
   
2. 2 #进入声子计算模块
   
3. 1 #选择slab模型
   
4. 1 #生成gamma点声子频率计算
   
5. How to calculate phonon frequency of adsorbed molecule at gamma:
   
6. * step 1: Perform a Self-Consistent Field ground-state calculation at the equilibrium structure using the pw.x program.
   
7.    Running cmd: mpirun -np 16 pw.x -i xxx.scf.in > xxx.scf.out
   
8. * step 2: calculate phonon frequency by ph.x code:
   
9.    Running cmd: mpirun -np 16 ph.x -i xxx.ph.in > xxx.ph.out
   
10. 
    
11. Input the atom index to be used in the linear response calculation, e.g. 3,5,10-15
    
12. 37-39 #这里仅仅对H2O进行线性响应计算
    
13. The input file has been generate at current folder!
    

复制代码

计算完成后，在TiO2_H2O.ph.out文件末尾就输出了H2O分子的3N个振动频率。为了计算吸附的H2O的热力学量，接下来我们继续启动QEtoolkit来制作Shermo的输入文件。

1. ./QEtoolkit TiO2.relax.in
   
2. 2 #ph.x模块
   
3. 1 #slab模型相关的输入文件生成
   
4. 2 #生成Shermo的输入文件以进行吸附分子热力学性质计算
   
5. TiO2_H2O.ph.out #还需要读取声子计算完后的输出文件。如果用的熟悉的话，直接./QEtoolkit TiO2_H2O.relax.in TiO2_H2O.ph.out 运行会更加方便。

复制代码

1. *E
   
2.   Input electronic energy in Hartree unit rather than Ry unit.
   
3. *wavenum
   
4. 172.812451
   
5. 268.561946
   
6. 288.998426
   
7. 369.409071
   
8. 564.328763
   
9. 662.637607
   
10. 1524.809370
    
11. 3699.021335
    
12. 3804.779306
    
13. *atoms
    
14. O       15.999   3.176546000    3.295033000    9.760404000
    
15. H       1.008    3.284608000    2.389300000    10.392720000
    
16. H       1.008    3.126065000    4.173219000    10.437422000
    
17. *elevel
    
18. 0.00    1
    

复制代码

运行Shermo的命令同样也被输出到了屏幕上，注意这里有个imode=1 的选项，是用来不计算吸附分子的平动和转动对热力学量的贡献的。

2.6 气相CO2分子的热力学性质计算

  首先建个边长为15埃的盒子，里面放个CO2，然后用pw.x进行自洽计算。接下来，如果你不清楚怎么计算气相分子的gamma点声子频率的话，那就可以进入相应模块的输入文件生成功能，屏幕上就会打印出计算的步骤。比如这里的CO2的频率及热力学量计算，运行./QEtoolkit CO2.scf.in，我们进入2 ph.x模块，然后选择2 气相分子，接着选1，输入文件就生成到了当前目录下，按照提示一步一步运行即可。

1. How to calculate phonon frequency of gaseous molecule at gamma:
   
2. * step 1: Perform a Self-Consistent Field ground-state calculation at the equilibrium structure using the pw.x program.
   
3.    Running cmd: mpirun -np 16 pw.x -i xxx.scf.in > xxx.scf.out
   
4. * step 2: calculate phonon frequency by ph.x code:
   
5.    Running cmd: mpirun -np 16 ph.x -i xxx.ph.in > xxx.ph.out
   
6. * step 3: Impose the Acoustic Sum Rule (ASR) by dynmat.x code:
   
7.    Running cmd: dynmat.x -i dynmat.in > dynmat.out

复制代码

对于气相分子， 我们必须将平动和转动对振动的贡献给投影掉，所以比起案例2.5中的吸附分子多出了在步骤2和3中多出了需要施加ASR的过程。第三个步骤计算完后，我们就会得到3N-5/6个频率，即对于线性分子是3N-5个频率，对于非线性分子是3N-6个频率，此例中对于CO2我们将会得到4个频率，可以在dynmat.mold中查看。
  对于CO2的热力学性质计算，我们同样是制作Shermo的输入文件，步骤如下所示：

1. ./QEtoolkit CO2.scf.in
   
2. 2 #进入ph.x模块
   
3. 2 #选择气相分子
   
4. 2 #制作Shermo的输入文件，此时将会自动读取dynmat.mold文件，该文件里面包含了CO2的频率信息
   

复制代码

当前目录下产生了CO2.shm文件，手动把里面第一行能量改为hartree单位制下的能量，然后运行Shermo CO2.shm -T 298.15 -P 1 -imode 0即可

1. *E
   
2. -38.02142
   
3. *wavenum
   
4.   870.04
   
5.   870.04
   
6. 1019.20
   
7. 1734.77
   
8. *atoms
   
9. C       12.011   6.000000000    6.000000000    6.000000000
   
10. O       15.999   7.258400000    6.000000000    6.000000000
    
11. O       15.999   4.741600000    6.000000000    6.000000000
    
12. *elevel
    
13. 0.00    1
    

复制代码

四、总结
    目前版本只是添加了一些最最常用功能，在以后的版本中会陆续添加比如SOC，RISM等等，MD那块我没有打算添加，做MD的话CP2K是最好的选择。计算原理这方面我还是一只菜菜鸟，仍有太多的知识需要学习，所以大家要是使用中遇到某些问题或者是我写错的地方，欢迎一起交流。

waitingseven

多谢开发和分享，我想问一下，QE里的SSSP赝势怎么样，后续会支持吗

丁越

waitingseven 发表于 2022-8-8 20:28
多谢开发和分享，我想问一下，QE里的SSSP赝势怎么样，后续会支持吗

这个我没用过唉，我都是用GBRV和PSlibrary中的赝势。你想用的话在脚本里自己定义一套SSSP赝势的数组就行了，然后照着里面的条件判断再改一改就可以用了。

waitingseven

丁越 发表于 2022-8-8 21:21
这个我没用过唉，我都是用GBRV和PSlibrary中的赝势。你想用的话在脚本里自己定义一套SSSP赝势的数组就行 ...

好的，多谢了。QE的赝势版本太多了，有时候不知道怎么选了，SSSP里每种元素都只有一个赝势，还有推荐的截断能，感觉会方便一些

丁越

waitingseven 发表于 2022-8-8 21:26
好的，多谢了。QE的赝势版本太多了，有时候不知道怎么选了，SSSP里每种元素都只有一个赝势，还有推荐的截 ...

工具包中17选项可以对截断能进行测试，每个体系计算之前都应该实际测试以下为好

flyingchow

丁老师您好，我刚刚用WSL运行了您的脚本，但是会报错，请问大概是什么原因呢？谢谢

QEtoolkit-2.sh: line 2: \r': command not found   
QEtoolkit-2.sh: line 11:\r': command not found
QEtoolkit-2.sh: line 22: syntax error near unexpected token `$'{\r''   
'Etoolkit-2.sh: line 22: `function main_menu (){   

丁越

flyingchow 发表于 2022-8-8 23:13
丁老师您好，我刚刚用WSL运行了您的脚本，但是会报错，请问大概是什么原因呢？谢谢

你是怎么运行的？比如你要使用0功能生成pw.x输入文件，就得载入如cif，gjf之类含有晶胞信息的文件，./QEtoolkit xxx.cif

flyingchow

按照这个来运行还是有报错，后来我shile “bash QEtoolkit-2.sh”，还是报错。

gog

Fedora 36报错同楼上
$ ~/QEtoolkit-2.sh
bash: /home/az/QEtoolkit-2.sh: /bin/bash^M: bad interpreter: No such file or directory。


ubuntu  22.04报错了

sh QEtoolkit-2.sh
QEtoolkit-2.sh: line 2: $'\r': command not found
QEtoolkit-2.sh: line 11: $'\r': command not found
>===============================================================================<
>-                                                                             -<
>-                                 QEtoolkit                                   -<
>-                                                                             -<
>-                              Author: Yue Ding                               -<
>-                                                                             -<
>-                              Version 2.0 (dev)                              -<
>-                                                                             -<
>-           Contact me if you have any suggestions: QQ 3305507756             -<
>===============================================================================<
QEtoolkit-2.sh: line 22: syntax error near unexpected token `$'{\r''
'Etoolkit-2.sh: line 22: `function main_menu (){

snljty2

flyingchow 发表于 2022-8-9 00:42
按照这个来运行还是有报错，后来我shile “bash QEtoolkit-2.sh”，还是报错。

1. dos2unix QEtoolkit-2.sh

复制代码

或者

1. cat QEtoolkit-2.sh | tr -d '\r' > QEtoolkit-2.sh.bak && mv QEtoolkit-2.sh.bak QEtoolkit-2.sh

复制代码

snljty2

gog 发表于 2022-8-9 07:23
Fedora 36报错同楼上
$ ~/QEtoolkit-2.sh
bash: /home/az/QEtoolkit-2.sh: /bin/bash^M: bad interprete ...

1. dos2unix QEtoolkit-2.sh

复制代码

或者

1. cat QEtoolkit-2.sh | tr -d '\r' > QEtoolkit-2.sh.bak && mv QEtoolkit-2.sh.bak QEtoolkit-2.sh

复制代码

gog

snljty2 发表于 2022-8-9 09:07
或者

it is all right, now. thks

tianmafei

gog 发表于 2022-8-9 09:33
it is all right, now. thks

请教楼上，这个问题是如何解决的?我也遇上了，不知道如何办。

gog

tianmafei 发表于 2022-9-2 09:37
请教楼上，这个问题是如何解决的?我也遇上了，不知道如何办。

见11楼

nianbin

可以联系卢老师整合进Multiwfn 这样就更好了 哈哈哈