3.3.1 结构优化参数介绍及实例#

(1)首先你需要一个POSCAR,内容如下:

Si8
1.0
5 0.0000000000 0.0000000000
0.0000000000 5 0.0000000000
0.0000000000 0.0000000000 5
Si
8
Direct
0.250000000 0.250000000 0.250000000
0.000000000 0.500000000 0.500000000
0.250000000 0.750000000 0.750000000
0.000000000 0.000000000 0.000000000
0.750000000 0.250000000 0.750000000
0.500000000 0.500000000 0.000000000
0.750000000 0.750000000 0.250000000
0.500000000 0.000000000 0.500000000

注意,这里的晶格常数是人为随意取的一个较为合理的值(为了演示用)。

(2)结构优化采用的INCAR内容如下:

ISTART = 0
ICHARG = 2
LCHARG = .FALSE.
LWAVE = .FALSE.
#PREC = Accurate
ENCUT = 400
ISPIN = 1
LREAL = A

# Ionic Relaxation
IBRION = 2
ISIF = 3
PSTRESS = 0
POTIM = 0.3
NSW = 100
EDIFFG = -0.01

# DOS Related values
ISMEAR = 0
SIGMA = 0.05
NEDOS = 2000
NELM = 60
EDIFF = 1E-06
#LORBIT = 11

NCORE = 4
KPAR = 4

相比于前面的静态自洽计算,结构优化多了IBRION,ISIF,NSW,POTIM,EDIFFG,PSTRESS等参数。

首先需要强调,这里对于Si来说,不需要考虑磁性、强关联(+U)、SOC、VDW等效应,也不采用高阶泛函(如PBEsol、SCAN等),因此INCAR中暂时不写出这些参数。其次我们对这些参数做一些说明:

PREC=Accurate,涉及到力的计算时(结构优化/声子频率等)取Accurate,其他计算可以注释掉。在这里我们注释掉了以减少计算量,可以自行比较注释前后的结果差异。

控制结构优化的主要参数:

NSW,最大离子步,当300个离子步还没有达到收敛标准就自动停止。

POTIM,控制结构优化过程中的原子移动距离,单位Å,默认结构优化时取0.5

IBRION,控制结构优化算法,在结构优化时可以取1 2 3三种,一般取2

ISIF,控制结构优化过程是否变晶格/胞。常用的是ISIF=2和3;2对应固定晶格优化原子位置(晶格长度不变),3是同时优化晶格与原子位置(晶格长度也会变化)。

EDIFFG= -0.01,单位是eV/Å,当每个原子在xyz方向的受力都小于该值时,优化停止。大部分情况下-0.01已经足够(对于一些金属表面吸附分子的问题,由于原子数很多,可以适当放宽精度,EDIFFG用-0.05);在算声子谱前的结构优化建议这个值设置为-0.001。

(EDIFFG可以取正值,那个对应的是前后两次优化的能量,这种情况用得不多)

请自行查阅官方手册的参数介绍:

https://www.vasp.at/wiki/index.php/Category:INCAR_tag

最后两行多了NCORE和KPAR参数,这两个参数是加速计算的参数,还有NPAR、KPAR参数。一般设置这两个参数比默认情况下,计算速度会更快。

https://mp.weixin.qq.com/s/uiWuvILE1WeiB5SC3_Frvg

(3)KPOINTS采用vaspkit-102生成,密度取0.025

K-Spacing Value to Generate K-Mesh: 0.025
0
Monkhorst-Pack
8 8 8
0.0 0.0 0.0

(4)POTCAR用vaspki-103生成。

然后确保这四个文件以及任务提交脚本在同一个路径下,通过脚本提交计算任务。

当任务顺利结束时,你应该要在result文件一行最后看到如下内容:

reached required accuracy - stopping structural energy minimisation

或者在OUTCAR最后看到如下信息:

reached required accuracy - stopping structural energy minimisation
    00P+    cpu time    0.2885: real time     0.3128
   4ORBIT:  cpu time    0.0000: real time     0.0000
total amount of memory used by VASP MPI-rank0 33717. kBytes
                            二二二二二二二二二二二二二二二二
  base            30000.kBytes
  nonlr-proj:       237. kBytes
  fftplans          583. kBytes
  grid             822. kBytes
  one-center:        6. kBytes
  wavefun          2069. kBytes
General timing and accounting informations for this job:
                Total CPU time used (sec):      2.893
                          User- time sec
                        System time (sec):      0.295
                       Elapsed time (sec):      3.947

接下来我们讲解一下结构优化计算过程

首先是result文件,会看到左边栏目有一系列的1F,2F等;每一个F代表一次循环,这里叫离子步,每一次离子步中的电子步都收敛了。这里一共有4F,即进行了四个离子步最终结构优化完成。

       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.453766044760E+02   -0.59780E+01   -0.69633E+01   982   0.357E+01    0.857E+00
DAV:   2    -0.440222202973E+02    0.13544E+01   -0.22006E+00  1010   0.690E+00    0.560E+00
DAV:   3    -0.431924393246E+02    0.82978E+00   -0.15979E+00   928   0.508E+00    0.165E+00
DAV:   4    -0.431944690917E+02   -0.20298E-02   -0.60683E-02  1040   0.121E+00    0.124E+00
DAV:   5    -0.431946622402E+02   -0.19315E-03   -0.11449E-03  1102   0.191E-01    0.987E-01
DAV:   6    -0.431952693995E+02   -0.60716E-03   -0.53724E-03   934   0.329E-01    0.330E-01
DAV:   7    -0.431952497998E+02    0.19600E-04   -0.91448E-05  1082   0.501E-02    0.233E-01
DAV:   8    -0.431952610584E+02   -0.11259E-04   -0.62956E-05   988   0.375E-02    0.206E-02
DAV:   9    -0.431952660620E+02   -0.50036E-05   -0.87567E-06   994   0.153E-02    0.404E-02
DAV:  10    -0.431952669714E+02   -0.90938E-06   -0.26879E-07   680   0.267E-03
   2 F= -.43195267E+02 E0= -.43195267E+02  d E =-.379667E+01
 trial-energy change:   -3.796667 1 .0rder     -4.755672 2-11.597769   2.086425
 step:   0.7740(harm=  0.8475) dis= 0.00000   next Energy=   -43.442139 (dE=-0.404E+01)
 bond charge predicted
       N       E                     dE             d eps       ncg     rms         rms(c)
DAV:   1     -0.434981277992E+02  -0.30286E+00   -0.39297E+00   982   0.901E+00    0.197E+00
DAV:   2    -0.434301857375E+02    0.67942E-01   -0.15027E-01   994   0.173E+00    0.120E+00
DAV:   3    -0.433890179907E+02    0.41168E-01   -0.71980E-02   910   0.110E+00    0.378E-01
DAV:   4    -0.433892200243E+02   -0.20203E-03   -0.22764E-03  1062   0.257E-01    0.326E-01
DAV:   5    -0.433891821565E+02    0.37868E-04   -0.24108E-04   990   0.908E-02    0.118E-01
DAV:   6    -0.433891678292E+02    0.14327E-04   -0.64030E-06  1068   0.164E-02    0.294E-02
DAV         -0.433891686606E+02   -0.83140E-06   -0.83507E-07   632   0.583E-03
         -.43389169E+02 E0= -.43389169E+02  d E =-.399057E+01
 curvature:  -0.33 expect dE=-0.183E-02 dE for cont 1inesearch -0.183E-02
 ZBRENT: interpolating
 opt :   0.7986 next Energy=    -43.392324 (dE=-0.399E+01)
 bond charge predicted
       N       E                     3P             d eps       ncg     rms          rms(c)
DAV:   1     -0.433931441495E+02  -0.39763E-02   -0.46061E-02   988   0.964E-01    0.222E-01
DAV:   2    -0.433922979272E+02    0.84622E-03   -0.15146E-03  1008   0.182E-01    0.137E-01
DAV:   3    -0.433917817837E+02    0.51614E-03   -0.88557E-04   898   0.122E-01    0.409E-02
DAV:   4    -0.433917824268E+02   -0.64308E-06   -0.27617E-05  1064   0.278E-02    0.347E-02
DAV:   5    -0.433917752248E+02    0.72020E-05   -0.20199E-06   730   0.854E-03    0.490E-03
DAV:        -0.433917747946E+02    0.43016E-06   -0.17689E-07   542   0.273E-03
   4 F=  .43391775E+02 E0= -.43391775E+02  d E =-.399317E+01
 curvature:  -0.56 expect dE=-0.213E-04 dE for cont 1inesearch -0.213E-04
 trial: gam= 0.00000 g(F)= 0.141E-59 g(S)= 0.380E-04 ort =-0.210E-01 (trialstep = 0.960E+00)
 search vector abs. value=  0.380E-04
 reached required accuracy - stopping structural energy minimisation

现在我们查看relax过程中每一次离子步的体系的压强。

用如下命令grep ' in kB ' OUTCAR 屏幕会输出OUTCAR中的in kB关键词(自行在屏幕输入命令内容,直接复制可能有格式问题)

in kB     418.96808    418.96808   418.96808     0.00000    -0.00000      0.00000
in kB     -53.89764    -53.89764   -53.89764     0.00000     0.00000     -0.00000
in kB       7.10339      7.10339     7.10339    -0.00000     0.00000     -0.00000
in kB      -0.57876     -0.57876    -0.57876    -0.00000     0.00000     -0.00000

这里输出了四行,对应四个离子步的压强,从左到右依次是Pxx,Pyy,Pzz,Pxy,Pyz,Pzx,可以自行翻看OUTCAR文件找到这些内容相应的位置。第一个离子步输出的体系压强是418kB,即42GPa(10kBar=1GPa),意味着当前晶格常数太小,体系会膨胀(有一个正压强)。最后压强变成了0(预设PSTRESS=0kBar),晶格常数到了5.47Å。

同时我们还关心每个离子步中每个原子的受力情况,

自行翻看OUTCAR,找到这部分内容:

POSITION                                       TOTAL-FORCE (eV/Angst)
     1.36807      1.36807      1.36807         0.000000     0.000000      -0.000000
     0.00000      2.73614      2.73614         0.000000     0.000000      0.000000
     1.36807      4.10421      4.10421         0.000000     -0.000000    -0.000000
     0.00000      0.00000     -0.00000        -0.000000     -0.000000     0.000000
     4.10421      1.36807      4.10421        -0.000000     0.000000      0.000000
     2.73614      2.73614      0.00000        -0.000000     0.000000      0.000000
     4.10421      4.10421      1.36807        -0.000000     -0.000000    -0.000000
     2.73614     -0.00000      2.73614        -0.000000     -0.000000     0.000000
   total drift:                               -0.000000     -0.000000     -0.000000

右边八行就是八个原子的受力,分别是xyz三个方向。这里每个原子受力都是0,是因为原子位置非常规整,没有偏离理想格点。也可以使用awk命令抓取力的部分:awk '/POSITION/,/total drift:/' OUTCAR (自行在屏幕输入命令内容,直接复制可能有格式问题)

结构优化的目的是为了得到合适的晶格常数和原子位置。该信息在CONTCAR文件中。我们查看CONTCAR文件得到relax后的结构的晶格常数,可以看到晶格常数变成了5.47Å。CONTCAR文件格式和POSCAR一样,因此可以直接把relax成功后输出的CONTCAR文件复制为POSCAR文件(新建文件夹操作,不要覆盖原来的relax中的POSCAR),用于后续其它计算。

Si8
   1.00000000000000
     5.4722791521771645   -0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.4722791521771645   -0.0000000000000000
     0.0000000000000000   -0.0000000000000000    5.4722791521771645
   Si
     8
Direct
  0.2500000000000000  0.2500000000000000  0.2500000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.2500000000000000  0.7500000000000000  0.7500000000000000
  0.0000000000000000  0.0000000000000000 -0.0000000000000000
  0.7500000000000000  0.2500000000000000  0.7500000000000000
  0.5000000000000000  0.5000000000000000  0.0000000000000000
  0.7500000000000000  0.7500000000000000  0.2500000000000000
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

XDATCAR文件里面有每一次离子步的结构信息,可以自行打开查看relax过程中的晶格常数变化。OVITO软件可以直接查看原子位置和轨迹,请自行网上下载(https://www.ovito.org/),并自学或向同门请教如何使用基本功能。