2.3.2 静态自洽计算的输出文件#

在这里,你需要跟着教程初步了解VASP的主要输出文件及其含义。

提交计算任务之后,如果任务正常开始计算时,我们会看到当前路径下多了很多文件,用ls命令查看

(base) [zpwu@mu003 Si]$  ls
2617356.err  DOSCAR     KPOINTS  POTCAR         vasprun.xml
2617356.out  EIGENVAL   OSZICAR  REPORT         WAVECAR
CHG          hostfile   OUTCAR   result         XDATCAR
CHGCAR       IBZKPT     PCDAT    SYMMETRY
CONTCAR      INCAR      POSCAR   vasp641.slurm

在这里,我们重点需要关注的是result、OUTCAR、OSZICAR这三个文件,除了第一个文件名字取决于脚本里的写法,其余文件名字都是VASP程序固定的命名。

(1)result文件

Result文件最上面几行会给出此次计算调用的计算资源(核数)、vasp版本等。之后会进入自洽迭代步骤(俗称,电子步)。请注意,result文件电子步前面部分可能有warning,一般可以不管,感兴趣可用自行去阅读这些warning内容。但是如果最后出现了warning,通常意味着计算可能有问题。这一点我们会在第三章详细说明。

(2)OUTCAR文件

VASP的全部计算信息都会在这个文件显示。输入cat OUTCAR并回车,然后借助鼠标滚轮上滑从第一行开始看。OUTCAR中前面依次会有vasp的版本信息,计算资源(核数),INCAR内容,POTCAR主要内容(这里不做展示)。

继续往下,看到了结构信息

  Lattice vectors:
 A1      5.4687280655,   0.0000000000,  0.0000000000)
 A2 =    0.0000000000,  5.4687280655    0.0000000000)
 A3 =    0.0000000000,  0.0000000000,   5.4687280655)
Analysis of symmetry for initial positions (statically):
 Subroutine PRIcEL returns following result:
  LATTYP: Found a face centered cubic cell.
 ALAT            5.4687280655
  Lattice vectors:
 A1 =    2.7343640327,  -2.7343640327,  0.0000000000)
 A2 =    2.7343640327,  0.0000000000,   2.7343640327)
 A3 =    0.0000000000,  -2.7343640327,  2.7343640327)
   4 primitive cells build up your supercell.

注意,这里vasp默认对结构寻找了对称性,因此找到了一个更小单元的胞进行计算,这个小胞里只有两个原子。

                                   Primitive cell
volume of cell :      40.8883
direct lattice vectors                    reciprocal lattice vectors
   2.734364033-2.734364033   0.000000000     0.182857876 -0.182857876 -0.182857876
   2.7343640330.000000000    2.734364033     0.1828578760.182857876   0.182857876
   0.000000000 -2.734364033  2.734364033    -0.182857876-0.182857876  0.182857876
length of vectors
   3.866974700 3.866974700   3.866974700     0.316719132 0.316719132  0.316719132
position of ions in fractional coordinates (direct lattice)
   0.750000000  0.750000000  0.750000000
   0.000000000  0.000000000  0.000000000

以及K点的信息

Subroutine IBzKPT returns following result:
Found      4 irreducible k-points:
Following reciprocal coordinates:
           Coordinates              Weight
 0.125000  0.125000 0.125000      8.000000
 0.375000  0.125000 0.125000     24.000000
 0.375000  0.375000 0.125000     24.000000
 0.375000  0.375000 0.375000      8.000000
Following cartesian coordinates:
           Coordinates              Weight
 0.022857  0.022857 0.022857      8.000000
 0.068572  0.022857 0.022857     24.000000
 0.068572  0.068572 0.022857     24.000000
 0.068572  0.068572 0.068572      8.000000

回想我们设置的KPOINTS文件,取的k点应该有64个。由于体系具有对称性,很多k点是等价的,实际只需要计算那些不可约的点,这里只有4个。在IBZKPT文件中也可以查看到这些不可约k点。当我们把每个k点的权重Weight相加,就得到了64个k点。在k点信息之后,还会给出完整的INCAR参数,因此实际计算中INCAR取的参数值也可以在OUTCAR这里查看。

SYSTEM = unknown system
POSCAR = Si8
Startparameter for this run:
  NWRITE =     2    write-flag & timer
  PREC  = normal   normal or accurate (medium, high low for compatibility)
  ISTART           job   :0-new  1-cont  2-samecut
  ICHARG       2    charge: 1-file 2-atom 10-const
  ISPIN        1    spin polarized calculation?
  LNONCOLLINEAR =    F non collinear calculations
  LSORBIT =     F    spin-orbit coupling
  INIWAV:      1    electr: 0-lowe 1-rand 2-diag
  LASPH        F    aspherical Exc in radial PAw
Electronic Relaxation 1
  ENCUT =  400.0 eV 29.40 Ry    5.42 a.u.  8.92 8.92  8.92*2*pi/ulx,y,z
  ENINI =  400.0     initial cutoff
  ENAUG =  322.1 eV augmentation charge cutoff
  NELM        60;  NELMIN= 2; NELMDL= -5     # of ELM steps

继续下滑,我们会看到这部分内容,左边就是4个不可约k点的坐标,最右边是每个k点对应的波函数所需要的平面波(可以理解为需要这么多平面波将该电子态的波函数做展开,这些平面波数量直接取决于INCAR中设置的ENCUT参数,参数越大,考虑的平面波越多,结果也越精确,代价就是贵(对于二维材料需要在z方向加真空层,真空层越大,所需要的平面波越多,后面我们会进一步了解)。

k-point  1     0.1250 0.1250 0.1250 plane waves: 2969
k-point  2     0.3750 0.1250 0.1250 plane waves: 2962
k-point  3:    0.3750 0.3750 0.1250 plane waves:  2970
k-point  4:    0.3750 0.3750 0.3750 plane waves: 2983
maximum and minimum number of plane-waves per node   2983    2962
maximum number of plane-waves:   2983
maximum index in each direction:
  IXMAX=   8   IYMAX=   8  IZMAX=    8
  IXMIN=  6-   IYMIN=  6-  IZMIN=   6-
serial  3DFFT for wavefunctions
parallel 3D FFT for charge:
   minimum data exchange during FFTs selected (reduces bandwidth)

继续下滑,我们可以看到:

                                  Iteration     1(  1)
POTLOK: cpu time      0.0075: real time    0.0130
SETDIJ: cpu time      0.0011: real time    0.0011
 EDDAV: cpu time      0.0504: real time    0.1203
   DOS: cpu time      0.0007: real time    0.0007
  LOOP: cpu time      0.0597: real time    0.1352

这是程序开始进行自洽迭代计算了,括号里的(1)意思是第一次循环,叫做第一个电子步。每次电子步算完,都会输出能量信息

Free energy of the ion-electron system (eV)
 alpha Z        PSCENC          13.17272099
 Ewald energy   TEWEN         -907.77962997
 -Hartree energ DENC   =       -45.29355852
 -exchange      EXHF   一         0.00000000
 -v(xc)+e(xc)   XCENC          -74.74068619
 PAW double counting           704.99039692     -564.24414232
 entropy T*S    EENTRO =        -0.00205010
 eigenvalues    EBANDS          73.94756885
 atomic energy  EATOM  =       824.52596983
 Solvation  Ediel_sol            0.00000000
 free energy    TOTEN           24.57658948 eV
 energy without entropy =       24.57863958  energy(sigma->0) :      24.57761453

经历多次电子步迭代,当前后两次能量差达到我们INCAR中设置的收敛阈值EDIFF时,程序停止迭代,开始计算电子结构和原子受力等信息。

比如这里输出了费米能级E-fermi和每个k点在每个态band(能带)上的电子占据occupation情况。

E-fermi :   5.6033     XC(G=0): -9.3267      alpha+bet :-11.7317
Fermi energy:         5.6032749100
k-point     1:        0.1250    0.1250   0.1250
 band No.  band energies     occupation
     1      -5.9979      2.00000
     2      -2.9678      2.00000
     3      -2.9678      2.00000
     4      -2.9678      2.00000
     5      -1.1620      2.00000
     6      -1.1620      2.00000
     7      -1.1620      2.00000
     8       2.2585      2.00000
     9       2.2585      2.00000

继续下滑,会看到计算达到收敛相关的信息输出

                    aborting loop because EDIFF is reached
CHARGE:  cpu time     0.0048: real time      0.0048
FORLOC:  cpu  time    0.0003: real time      0.0003
FORNL    cpu time     0.0117: real time      0.0117
STRESS:  cpu time     0.0319: real time      0.0320
FORCOR:  cpu time     0.0029: real time      0.0029
FORHAR:  cpu time     0.0011: real time      0.0011
MIXING:  cpu time     0.0003: real time      0.0003
OFIELD:  cpu time     0.0000: real time      0.0000

随后开始计算压强和受力

FORcE on cell =-STRESS in cart. coord.  units (eV):
Direction    XX                      ZZ         XY          YZ          ZX
Alpha Z    13.17272    13.17272    13.17272
Ewald    -302.59373  -302.59373  -302.59373    -0.00000    -0.00000    -0.00000
Hartree    20.23882    20.23882    20.23882    -0.00000    -0.00000    -0.00000
E(xc)    -100.98605  -100.98605  -100.98605    -0.00000    -0.00000    0.00000
Local    -118.08797  -118.08797  -118.08797     0.00000    -0.00000    0.00000
n-local   309.70644   312.31673   307.06329    -3.39890    -2.83265    -2.51354
augment   -46.57837   -46.57837   -46.57837     0.00000     0.00000    0.00000
Kinetic   218.93795   228.79224   227.35762    -4.58272    -6.30846    -4.07425
Fock        0.00000     0.00000     0.00000     0.00000     0.00000    0.00000
Total      -0.10981    -0.10981    -0.10981     0.00000     0.00000    -0.00000
in kB      -1.07575    -1.07575    -1.07575     0.00000     0.00000    -0.00000
external pressure         -1.08 kB  Pullay stress =       0.00 kB

其中,Total后面的六个数字是位力(Virial),单位是eV,六个数值分别对应XX, YY, ZZ, XY, YZ, ZX六个分量。请自行补“位力”相关的物理背景知识。in kB后面的六个数字也是Virial,单位是 kbar(1GPa=10kB,1B就是1bar,即一个标准大气压),数值上直接等于Total那一栏的数字除以当前这一帧结构的体积,注意单位换算。这个单位下的位力,可以简单理解为当前结构的静态压强。可以发现此时整个体系会有个负压强external pressure = -1kB(对应0.1GPa,1GPa=10kB,1B就是1bar,即一个标准大气压),说明晶格在此时有自发向内收缩的趋势;如果这个数值是正的,意味着晶格在当前压强下有膨胀趋势。对于晶体来说,0.1GPa对晶格常数、带隙等性质影响非常弱。

压强之后会输出原子受力信息:

POSITION                                       TOTAL-FORCE (eV/Angst)
     1.36718      1.36718      1.36718         0.000000      000000*0     -0.000000
     0.00000      2.73436      2.73436        -0.000000      0.000000      0.000000
     1.36718      4.10155      4.10155         000000*0     -0.000000     -0.000000
     0.00000      0.00000      0.00000         0.000000     -0.000000      0.000000
     4.10155      1.36718      4.10155        -0.000000      0.000000     -0.000000
     2.73436      2.73436      0.00000        -0.000000      000000*0      0.000000
     4.10155      4.10155      1.36718         0.000000     -0.000000     -0.000000
     2.73436      0.00000      2.73436        -0.000000     -0.000000      0.000000
   total drift:                               -0.000000     -0.000000      0.000000

由于我们的结构POSCAR中原子位置非常高对称,因此每个原子受力TOTAL-FORCE都为0。这只是我们当前例子比较特殊导致的巧合。

最后是关键的能量信息:

FREE ENERGIE OF THE ION-ELECTRON SYSTEM (ev)
free  energy   TOTEN          -43.39328400eV
energy  without entropy=      -43.39328400  energy(sigma->0)        43.39328400

这之后就是整个计算所运行的时间(1.78秒)

General timing and accounting informations for this job:
                Total CPU time used (sec):       1.783
                          User time (sec):       1.399
                        System time (sec):       0.384
                       Elapsed time (sec):       2.802
                 Maximum memory used (kb):     221260.
                 Average memory used (kb) :        N/A
                        Minor page faults:       21065
                        Major page faults:
               Voluntary context switches:        4565

当你看到OUTCAR最后出现这个信息时,意味着计算成功运行且结束了,但并不意味着自洽计算部分收敛了。需要结合前面提到的“EDIFF is reached”信息以及OSZICAR或者result文件中的电子步。

(3)OSZICAR,输入cat OSZICAR并回车

(base) [zpwu@muoo3 Si]$ cat OSZICAR
       N       E                     dE             d eps       ncg      rms          rms(c)
DAV:   1     0.245765894816E+02    0.24577E+02    -0.11363E+04  256    0.658E+02
DAV:   2    -0.432737377084E+02   -0.67850E+02    -0.67296E+02   384   0.786E+01
DAV:   3    -0.440666097061E+02    -0.79287E+00   -0.79287E+00   384   0.112E+01
DAV:   4    -0.440672202112E+02    -0.61051E-03   -0.61051E-03   512   0.371E-01
DAV:   5    -0.440672202424E+02    -0.31149E-07   -0.31149E-07   320   0.312E-03    0.948E+00
DAV:   6    -0.436028846451E+02    0.46434E+00    -0.30374E-01  416    0.144E+00    0.580E+00
DAV:   7    -0.433887736873E+02    0.21411E+00    -0.59457E-01   384   0.215E+00    0.326E-01
DAV:   8    -0.433920258292E+02    -0.32521E-02   -0.13379E-02   320   0.432E-01    0.155E-01
DAV:   9    -0.433931954126E+02    -0.11696E-02   -0.13474E-03  480    0.138E-01    0.106E-01
DAV:  10    -0.433932757836E+02   -0.80371E-04    -0.66202E-05   320   0.313E-02    0.792E-03
DAV:  11    -0.433932822670E+02    -0.64834E-05   -0.83780E-06  480    0.103E-02    0.568E-03
DAV:  12    -0.433932839156E+02    -0.16486E-05   -0.26691E-06   320   0.646E-03    0.227E-03
DAV:  13    -0.433932839974E+02    -0.81813E-07   -0.83379E-09  320    0.422E-04
   1 F=-.43393284E+02  E0= -.43393284E+02  dE  =-.200187E-14

说明程序运行了13个电子步达到收敛,每个电子步有一个dE,当dE小于INCAR中设置的EDIFF=1E-6时,scf结束。

并且输出了基态总能量E0,以及F还有一个d E=-0.829E-3(这个d E和另外一个dE不同,前者是以正值EDIFFG为收敛阈值的)

我们再查看一下OUTCAR中最后一个电子步给出的能量信息

Free energy of the ion-electron system (eV)
 alpha Z        PSCENC          13.17272099
 Ewald energy   TEWEN         -907.77962997
 -Hartree energ DENC           -60.71687757
 -exchange      EXHF             0.00000000
 -v(xc)+e(xc)   XCENC          -72.23214521
 PAW double counting          3311.55222047    -3171.62851930
 entropy T*S    EENTRO          -0.00000000
 eigenvalues    EBANDS          19.71297677
 atomic energy  EATOM          824.52596983
 Solvation  Ediel_sol            0.00000000
 free energy    TOTEN          -43.39328400eV
 energy without entropy =      -43.39328400  energy(sigma->0) =      -43.39328400

关注最下面的三个能量:free energy,energy without entropy和energy(sigma-0)。这里的free energy就是OSZICAR中的F,energy(sigma-0)就是OSZICAR的E0. Free energy是对蓝色框中的全部能量求和,其中EATOM来源于POTCAR(可以自行查看POTCAR前面几行)。E0就是在F的基础上扣掉了电子熵对能量的贡献(entropy TS)并且对展宽sigma取极限得到的,所以基态能量取E0。由于金刚石Si在常压下是半导体,因此不需要考虑电子熵的贡献问题(费米面上没有电子态占据);但如果是金属,电子熵影响就需要考虑了。无论如何,在常规的计算中,以E0即energy(sigma-0)作为体系的能量。特殊地,某些高温下的金属体系模拟,是需要考虑电子熵的。关于能量E0的理解可以阅读:https://mp.weixin.qq.com/s/1rsCjl-f_w6Es1wZ6kw_Dw

需要强调的是,这里直接算出来的能量可以简单理解为是电子和离子相互作用的势能,整个体系被认为是完全静止的。没有考虑原子振动对能量的贡献,这部分能量需要通过计算声子获得,超出了这一节范畴。

其余输出文件除了WAVECAR之后,其余都可以直接以文本形式打开查看。在大部分计算中,用到较多的输出文件除了上面三个result,OUTCAR,OSZICAR之外,还有CONTCAR,CHGCAR,WAVECAR(不能直接打开查看),DOSCAR,EIGENVAL,IBZKPT,XDATCAR(结构驰豫或者AIMD过程中的原子轨迹)。这些文件的内容请查阅官网,直接在vaspwiki搜索栏输入文件名(大写)即可。