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