使用 QE 进行能量自洽计算

pw.x 处理的计算包括以下7种类型，在输入文件中用 calculation 设置：

'scf'：自洽计算，self-consistent field，通过迭代的方式数值求解微分-积分方程（Kohn-Sham方程），迭代收敛以电荷的变化足够小为准，最终得到自洽电荷。

'nscf'：非自洽计算，nscf 计算常在k空间的网格上进行，网格要足够密以完成 k 空间上的积分，在DOS等计算需要更密的 k 点，这时在自洽电荷基础上，计算这些更多的 k 点，nscf计算保持自洽电荷不变。

'bands'：也是一种 nscf 计算，k 点按照三维 k 空间中的特殊路径选取。

'relax'：一系列 scf 计算，通过 Hellman-Feynman 力计算离子坐标驰豫（通过优化算法找到受力为零的结构），relax计算时固定cell不变。

'vc-relax': 允许cell变化的relax，通过应力的计算改变cell。

'md'：分子动力学，将电子对离子的作用看成离子感受到的势，根据势能和离子初始速度求解离子运动的经典力学方程。

'vc-md'：允许cell改变的md。

输入文件

在 QE 的输入文件中，有 NAMELISTS 和 INPUT_CARDS。

PWscf 中有三个强制 NAMELISTS：

1. &CONTROL：指定计算流，
2. &SYSTEM：指定系统，
3. &ELECTRONS：指定用于求解 Kohn-Sham 方程的算法。

还有另外两个 NAMELISTS：&IONS 和 &CELLS，需要根据计算指定。

PWscf 中的三个 INPUT_CARDS：ATOMIC_SPECIES、ATOMIC_POSITIONS 和 K_POINTS 是必需的。 在某些计算中还必须提供其他一些信息。

示例：单晶硅的自洽计算 pw.scf.silicon.in

&CONTROL
calculation='scf',   # self-consistent field 进行自洽计算
prefix='silicon',   # 输出文件的前缀
pseudo_dir='./pseudo/',   # 存放赝势的文件夹
outdir='./out/',   # 指向输出文件的文件夹
/
&SYSTEM
ibrav=2,   # 布拉维晶格类型 FCC结构
celldm(1)=10.2625,   # 晶格常数 （单位为波尔，1波尔=0.529A）
nat=2,   # 结构内原子个数
ntyp=1,   # 结构内原子种类数
ecutwfc=60.0,   # 波函数截断能 （单位为里德伯格，1Ry = 13.606 eV）
ecutrho=720.0,   # 电荷密度截断能
! occupations = 'smearing', smearing = 'gauss', degauss = 1.0d-9 # 高斯展宽
/
&ELECTRONS
mixing_beta=0.7,   # 自洽计算过程中前后电荷密度混合比，默认0.7
conv_thr=1d-8,   # 自洽收敛阙值，默认为1d-6
/
ATOMIC_SPECIES
SI 28.0855 Si.pbe-rrkj.UPF   # 赝势
ATOMIC_POSITIONS (alat)   # 原子位置（alat 表示采用分子坐标）
Si 0.00 0.00 0.00
Si 0.25 0.25 0.25
K_POINTS automatic  # K点自动分布
4 4 4 1 1 1

• ibrav - 取值 1-14 对应 14 个布拉维晶格，0对应自由晶格
• ecutwfc - 波函数截断能，选择不同的赝势，会有对应的建议值。
• ecutrho - 电荷密度截断能，与 ecutwfc 大小相关
• 赝势为 NC ，该值为 ecutwfc 的 4 倍
• 赝势为 US，该值为 ecutwfc 的 8-12 倍
• 赝势为 PAW，该值为 ecutwfc 的 4 倍
• mixing - 自洽计算过程中前后电荷密度混合比，默认0.7
• conv_thr - 自洽收敛阙值，默认1d-6
• 文件保存为 pw.scf.silicon.in
• 建议设置 occupations='smearing',smearing='gaussian'，使用高斯展宽。高斯展宽对导体和绝缘体是通用的，对于不同的材料调整degauss的值：（1）对于有带隙的材料，degauss=1d-9，也就是接近0就可以；（2）对于没有带隙的导体材料，degauss 从小尝试到大。

必须阅读 PWscf user manual 才能深入理解自洽计算的输入文件。

https://www.quantum-espresso.org/Doc/INPUT_PW.html

运行

mpirun -np 4 pw.x -inp pw.scf.silicon.in > pw.scf.silicon.out

输出文件

现在让我们看看输出文件 pw.scf.silicon.out 并看看如何达到收敛：

grep -e 'total energy' -e estimate pw.scf.silicon.out

之后我们会看到：

Self-consistent Calculation

iteration #  1     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  1.00E-02,  avg # of iterations =  2.0
Threshold (ethr) on eigenvalues was too large:
Diagonalizing with lowered threshold
Davidson diagonalization with overlap
ethr =  6.93E-04,  avg # of iterations =  1.0
total cpu time spent up to now is        0.3 secs
total energy              =     -15.83539933 Ry
estimated scf accuracy    <       0.06071141 Ry

iteration #  2     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  7.59E-04,  avg # of iterations =  1.0
total cpu time spent up to now is        0.4 secs
total energy              =     -15.83851631 Ry
estimated scf accuracy    <       0.00218630 Ry

iteration #  3     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  2.73E-05,  avg # of iterations =  2.3
total cpu time spent up to now is        0.4 secs
total energy              =     -15.83897965 Ry
estimated scf accuracy    <       0.00007075 Ry

iteration #  4     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  8.84E-07,  avg # of iterations =  2.6
total cpu time spent up to now is        0.4 secs
total energy              =     -15.83900178 Ry
estimated scf accuracy    <       0.00000211 Ry

iteration #  5     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  2.64E-08,  avg # of iterations =  3.6
total cpu time spent up to now is        0.4 secs
total energy              =     -15.83900299 Ry
estimated scf accuracy    <       0.00000012 Ry

iteration #  6     ecut=    16.00 Ry     beta= 0.70
Davidson diagonalization with overlap
ethr =  1.46E-09,  avg # of iterations =  2.7
total cpu time spent up to now is        0.4 secs

End of self-consistent calculation

值得注意的是，DFT 总能量的绝对值并不相对于真空参考，而是取决于所选的赝势。

!    total energy              =     -15.83900302 Ry
     estimated scf accuracy    <          3.2E-10 Ry

     The total energy is the sum of the following terms:
     one-electron contribution =       4.79863652 Ry
     hartree contribution      =       1.07565897 Ry
     xc contribution           =      -4.81353993 Ry
     ewald contribution        =     -16.89975858 Ry

     convergence has been achieved in   6 iterations

Tips

• 减少 mixing_beta 值，特别是当收敛能量周围存在振荡时。
• 增加能量和电荷密度截止值（确保它们足够）。
• 某些伪势文件有问题，可以尝试使用不同库中的伪势。
• conv_thr 的建议值：
  • 能量和特征值（scf 计算）1.0d-7 Ry
  • 受力（relax 计算）1.0e-8 Ry
  • 应力（vc-relax 计算）1.0d-9 Ry
• 能量和特征值（scf 计算）1.0d-7

输出文件上还打印了其他几个重要信息。

计算中使用的交换相关性：

Exchange-correlation= SLA  PZ   NOGX NOGC

其中 SLA — Slater交换； PZ — LDA 的 Perdew-Zunger 参数化； NOGX 和 NOGC 表示不考虑密度梯度。

计算中使用的平面波总数 (1067)：

Parallelization info
--------------------
sticks:   dense  smooth     PW     G-vecs:    dense   smooth      PW
Min         108     108     34                 1489     1489     266
Max         109     109     35                 1492     1492     267
Sum         433     433    139                 5961     5961    1067

Kohn-Sham 态的数量：

number of electrons       =         8.00
number of Kohn-Sham states=            8