变胞弛豫vc-relax和弛豫relax

@https://yyyu200.github.io/DFTbook/blogs/2019/04/01/HandsOn/#121-%E8%87%AA%E6%B4%BD%E8%AE%A1%E7%AE%97mathbfk%E7%82%B9%E7%9A%84%E8%AE%BE%E7%BD%AE

pw.x处理的计算包括的类型：

calculation=’scf’

自洽计算，通过迭代的方法数值求解KS方程，迭代收敛以电荷的变化足够小为标准，最终得到自洽电荷

calculation=’nscf’

非自洽计算，scf计算通常在k-mesh上，网格需要足够密来完成k空间上的积分，非自洽需要在自洽基础上取更密的k-mesh

calculation=’bands’

特殊的非自洽计算，k点常选择三维k空间特殊的路径

calculation=’relax’ 固定cell不变，只改变坐标位置

一系列scf计算，通过Hellman-Feynman力计算离子坐标弛豫（通过优化算法找到受力为零的结构）

calculation=’vc-relax’

变胞计算，通过应力的计算变胞

calculation=’md’

分子动力学，将电子对离子的作用看作是离子感受到的势场，根据势能和离子的初始速度求解离子运动的牛顿方程

calculation=’vc-md’

变胞的分子动力学

结构弛豫的输入文件

&CONTROL
    calculation='vc-relax', disk_io='low', prefix='pwscf',
    pseudo_dir='./', outdir='./tmp', verbosity='high'
    tprnfor=.true., tstress=.true., forc_conv_thr=1.0d-5
/
&SYSTEM
    ibrav= 0,
    celldm(1) = 1.8897261328856432, ! a.u. to Angst
    nat= 2, ntyp= 2,
    occupations = 'smearing', smearing = 'gauss', degauss = 1.0d-9
    ecutwfc= 50, ecutrho = 500,
/
&ELECTRONS
    electron_maxstep = 100
    conv_thr = 1.0d-9
    mixing_mode = 'plain'
    mixing_beta = 0.8d0
    diagonalization = 'david'
/
&IONS
    ion_dynamics='bfgs'
/
&CELL
    press_conv_thr=0.1
/
ATOMIC_SPECIES
  Si 28.08550 Si.UPF
  C  12.01070 C.UPF
CELL_PARAMETERS (alat=  1.88972613)
   2.174000000   2.174000000   0.000000000
   0.000000000   2.174000000   2.174000000
   2.174000000   0.000000000   2.174000000
ATOMIC_POSITIONS (crystal)
Si       0.000000000   0.000000000   0.000000000
C        0.250000000   0.250000000   0.250000000
K_POINTS {automatic}
  8 8 8 0 0 0

注意事项小结：

1.收敛标准需要足够提高

如forc_conv_thr,press_conv_thr, conv_thr在计算例子里常用的标准：

forc_conv_thr=1.0d-7
press_conv_thr=0.1
conv_thr=1.0d-10
 

2.尽可能提高截断能：

qe比较吃截断能，目前已知的一些材料，如FGT的ecutwfc= 120，Ni的ecutwfc= 50，总之，在计算资源充足的情况下要尽可能地提高截断能进行测试

3.smearing要分类讨论：

当材料是绝缘体，常用fixed或者gauss(此时degauss尽可能小一些，如1.0d-9)；当材料是导体或金属，常用mv，mp，但更推荐gauss！

4.degauss的测试：

针对3的讨论，无论什么才可以采用gauss，但需要测试degauss和k-mesh的匹配程度，在k-mesh取足够大(最小0.04)的情况下，测试degauss从1d-9~0.1(如常用的0.0001，0.0005，0.001，0.005，0.01，0.02，0.05，0.1等)

判断degauss合适的标准：

1)总能几乎不变,如 grep "!" pw.FGT.out

2) 输出中的“smearing contribution TS”足够小，建议小于1×10−4 Ry/atom

虽然degauss这个变量有时称为electronic temperature，但是，在scf计算过程中，设置展宽的目的并不代表温度，费米-狄拉克分布虽然是自然的选择，但是其取值的长尾效应，需要较大的degauss，设置smearing=’fermi-dirac’会造成收敛相对较慢，一般不推荐这种设置

k点和展宽的设置，不仅对于总能，而且对于所有需要在整个布里渊区积分的物理量的计算收敛有影响，所以，k点和展宽的选取，以所关心的物理量计算收敛为最终原则，同时考虑可接受的计算量

5.nbnd:无论增加k-mesh还是degauss，体系还是无法收敛，尝试增加nbnd

6.计算结束后，运行：

awk  '/Begin final coordinates/,/End final coordinates/{print $0}' vc.out

得到：

Begin final coordinates
     new unit-cell volume =    141.54631 a.u.^3 (    20.97500 Ang^3 )
     density =      3.17432 g/cm^3

CELL_PARAMETERS (angstrom)
   2.188890250   2.188890250   0.000000000
   0.000000000   2.188890250   2.188890250
   2.188890250  -0.000000000   2.188890250

ATOMIC_POSITIONS (crystal)
Si       0.000000000  -0.000000000   0.000000000
C        0.250000000   0.250000000   0.250000000
End final coordinates

即为弛豫后的晶格结构
vc-relax包括多个自洽的计算，查看自洽得到的总能量用：
grep '!' vc.out输出：

!    total energy              =     -22.67776498 Ry
!    total energy              =     -22.67819319 Ry
!    total energy              =     -22.67819977 Ry
!    total energy              =     -22.67819977 Ry
!    total energy              =     -22.67821030 Ry

从上可见5个自洽计算步骤后达到收敛，最后一个数值是结构弛豫计算得到的体系总能。
其中，第5步是在第4步压强收敛之后，用第4步的结构重新做了一次scf计算，由于CELL的变化，平面波的个数随之变化，最后第5步计算的压强和第4步会有一定的差异，如果差异较大，则说明截断能未收敛，需要增加截断能，或用最后的结构重新启动vc-relax