1.2 结构优化2

发表于 更新于 分类于

结构优化完成后，我们已经可以从输出文件得到很多性质。这里我们先关注能量。

VASP输出的能量是什么

在优化完的文件夹下输入

{.line-numbers}
1
2
[scf1355@ln151%nc-l 2h]$ tail -1 OSZICAR 
3 F= -.44829871E+02 E0= -.44829871E+02  d E =-.169526E-03  mag=    -0.0000

这里的F和下面free energy含义相同，E0代表体系内能。
另外在OUTCAR中还可以找到这些
{.line-numbers}
1
2
3
free  energy   TOTEN  =       -44.82987077 eV
energy  without entropy=      -44.82987077  energy(sigma->0) =      -44.82987077
atomic energy  EATOM  =      1542.31592482

free energy指的是亥姆霍兹自由能(Helmholtz free energy)，在VASP手册中被称为自由电子能(free electronic energy)，它可以用F表示： F = E-TS+PV .而F和E对应上面的F和E0。
TOTEN和without entropy分别代表考虑和未考虑熵的贡献时自由能的值。而energy(sigma->0)为展宽趋近于零时的自由能，可以外推得到： energy(sigma->0) = E0 =1/2（F+E）

其中上面的展宽对应于INCAR中设置的ISMEAR参数，详情可以查看手册或者参考 https://zhuanlan.zhihu.com/p/384626631。
https://www.bilibili.com/video/BV19L4y1i7mC/?vd_source=17084cc867ff3bb86398ff1a79b443f2
简单理解费米面内电子全占据，费米面外没有电子。这种阶梯函数就使得，随着K点数目变化，K点处积分值变化很大，不易收敛。于是我们需要一个连续函数，代替该阶梯函数进行积分，这种替代方法就叫smear。在VASP计算中，只有在ISMEAR=-5时，才能严格满足费米能级的定义；而在其他情况下，允许电子在费米能级附近的宽度内“波动”，即将阶梯函数替代为一个平滑的函数，从而在不破坏和的精度的情况下获得更快的收敛速度。

明白非常重要的一点，从OUTCAR中直接提取的能量在不同赝势下不具备重复性，是不能够直接写到论文中去的。一般计算相对能量的时候只要保持取一样的值即可 比较绝对能量没啥意义。

形成能

形成能很简单，参考高中化学形成焓的计算可以知道对于反应A+B=C+D，那么形成能=E(C)+E(D)-E(A)-E(B)
同样对于MoS2的生成焓我们可以知道： $\Delta$ H=E(MoS2)-E(Mo)-2E(S)
因此还需要计算Mo单质和S单质的单原子能量。计算过程和1.1完全一致，只需要获得结构，利用vaspkit生成KPOINTS和POTCAR文件，INCAR可以不用改。能量用F或者E0或者上面说的哪一项都可以，但要保持一致。代入计算即可。

除此之外还可以计算产生一个Mo的空位的形成能： $\Delta$ H=E(空位，MoS2)+E(Mo)-E（MoS2）
但是此时没有考虑缺陷的电荷。
可以计算分解能，剥离能，滑移能，表面能，偏聚能，吸附能等等。但是具体结构的优化会有一些处理上的不同。

单点能SP

参考： https://www.bilibili.com/read/cv10064345

单点计算是指对结果不做优化，直接输出能量。因此优化只会进行一步。

设置NSW = 0 也可以不设置，因为默认值就是0。 其他文件和参数不需要改变。提交即可。

对于一个结构在结构优化（NSW不等于0）时，第一步输出的1F和单点能的结果是一样的。

自洽计算

自洽计算其实也是单点计算，静态自洽计算前需要提供 较稳定体系 的晶格结构信息（即结构优化完的CONTCAR），从而通过电子自洽计算，通过自洽迭代求解薛定谔方程（微观中描述电子的状态，相当于宏观的运动方程）） 完整地计算出体系基态下费米能级（准确值，The fermi level location is accurate only in self-consistent calculation.）、 电子的波函数、电荷密度等信息，可以直接分析原子间的键合作用，也可以在非自洽之后进一步分析晶体的电子结构和材料的相关性质。

因此对于电子结构计算之前必须进行一步自洽计算优化电子分布，以便加速后续计算。

{.line-numbers}
1
2
3
NSW=0  不再进行原子迟豫
LWAVE=.TRUE  （需要在WAVECAR中输出波函数）
LCHARGE=.TRUE    （需要在CHG和CHGCAR中输出电荷密度）

气体优化

在涉及单质结构优化时，可能会涉及到气体分子的优化。

  • 一般情况下选择建一个相对大的格子，里面放入一个气体分子，格子的大小保证分子间不会影响。其他的和前面结构优化保持一致即可。
  • 另一个办法就是找数据库， ChemScpider。网址：http://www.chemspider.com

IVDW参数

有两大类的范德华斯作用修正：

  1. 基于半经验的，包括了D2, D3, D3-BJ, TS, TS+SCS等等，这些都是在常用的交换关联泛函比如PBE的基础上，考虑色散力的作用，在PBE计算出来的总能基础上增加了额外半经验项，这一项需要设置一些参数，至于是哪种半经验公式，由IVDW的设置来决定，具体的参考VASP的手册。
    范德华力计算方法，在DFT能量计算基础上增加范德华力修正。IVDW=10：DFT-D2方法；IVDW=11：DFT-D3方法。推荐首选更新的DFT-D3方法。
  2. 另一类就是vdW-DF，也就是修改交换关联项来考虑范德华修正。包括你所提到的optPBE-vdW, optB88-vdW, and optB86b-vdW等。它们的设置由修改GGA赋值，以及额外相关的一些参数来设置。 这些在vasp手册上有明确的说明。
    乱谈DFT-D：http://sobereva.com/83
    DFT-D色散校正的使用：http://sobereva.com/210
    大体系弱相互作用计算的解决之道：http://sobereva.com/214
    LUSE_VDW：https://blog.csdn.net/Yaris_liu/article/details/124683272 https://zhuanlan.zhihu.com/p/47619708
    包含vdW-DF，optPBE-vdW, optB88-vdW, optB88b-vdW, vdw-DF2。需要在所在的文件夹下放置一个名叫vdw_kernel.bindat的文件。
  • vdW-DF：采用该范德瓦尔斯修正需在INCAR中加上以下开关：
    {.line-numbers}
    1
    2
    3
    GGA=RE
    LUSE_VDW=.TRUE.
    AGGAC=0.0000
  • optB88-vdW：采用该范德瓦尔斯修正需在INCAR中加上以下开关：
{.line-numbers} GGA
1
2
3
4
PARAM1=0.1833333333
PARAM2=0.2200000000
LUSE_VDW=.TRUE.
AGGAC=0.0000
  • optB86b-vdW：采用该范德瓦尔斯修正需在INCAR中加上以下开关：
{.line-numbers} GGA
1
2
3
4
PARAM1=0.1234
PARAM2=1.0000
LUSE_VDW=.TRUE.
AGGAC=0.0000
  • vdW-DF2：采用该范德瓦尔斯修正需在INCAR中加上以下开关：
{.line-numbers} GGA
1
2
3
LUSE_VDW=.TRUE.
Zab_vdW=-1.8867
AGGAC=0.0000

在考虑DFT-D2或者DFT-D3计算的过程中可同时考虑自旋轨道耦合（SOC），但是在vdW-DF、optB88-vdW 、optB86b-vdW和vdW-DF2中无法同时考虑自旋轨道耦合进行计算。

DFT+U

  1. 为什么加U?
    对于含有d、f轨道电子的强关联体系(Hubbard模型)，电子间存在强烈的在位库仑相互作用，而交换相关泛函中的局域密度近似（LDA）或广义梯度近似（GGA）对电子之间的强在位库仑相互作用描述不准， 此时可将研究体系的交换相关泛函分成两部分计算，也就是用DFT+U的方法进行求解，公式如下：
    其中一部分电子用DFT算法（如LDA，GGA）等可以比较准确地描述；另外的d、f轨道电子通过引入Hubbard项得到正确的描述（U值就是考虑了同一个原子自旋相反的局域电子之间的库仑排斥）
  2. 如何加U?
    DFT计算时输入文件的U值由以下参数确定：
    {.line-numbers}
    1
    2
    3
    4
    5
    LDAU= .TRUE.|.FALSE.：开启/关闭+U功能，默认值为.FALSE.；
    LDAUTYPE=1|2|4指的是+U的类型，默认值是2；其中1为Liechtenstein等提出的旋转不变LSDA+U方法；2为 Dudarev等提出的简化 LSDA+U方法；4与1类似, 但不考虑LSDA交换劈裂；
    LDAUL=-1|1|2|3分别对应不加U、p、d、f轨道加U；
    LDAUU、LDAUJ分别指定电子库伦相互作用项和交换相互作用项（U和J值）；
    LMAXMIX =2/4/6：默认为2，加U计算时该值需大于轨道量子数，对于含有d轨道或f轨道电子的体系需对应增加至4或6。
    注意：
    1）必须为每一类原子指定一个U值，其顺序与结构文件POSCAR中元素排序一致；
    2）同一种元素在不同的体系中U值不同；
    3）U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上；当LDAUTYPE =2时，U-J的差值才有意义，即有效的U参数。
    该参数需要测试或者参考文献。

磁性体系

很多体系都有磁性，磁性体系的计算需要打开自旋。以FeO为例：

{.line-numbers}
1
2
ISPIN=2
MAGMOM=原子数目*初始磁矩，eg 3*2 4*5 （注：*号前后不要有空格，且对应顺序）

对于一些简单的磁性体系，我们可以直接使用ISPIN=2, MAGMOM不必进行设置。对于复杂的体系则需要设置。初始值不需要与实验值完全一致，一般设置1.5倍即可。

  1. 有4个Fe原子，如果考虑铁磁态，则可以设置 MAGMON= 4*3 4*0 （Fe初设的磁矩为3，O元素为0）

  2. 有4个Fe原子，如果考虑反铁磁态，则可以设置 MAGMON=3 -3 3 -3 0 0 0 0或者MAGMON=3 3 -3 -3 0 0 0 0

表面优化

1. Bulk计算，获得稳定构型。

注意LREAL、PREC、ENCUT、EDIFF、ISMEAR保持一致。

2. 考虑不同晶面、层数。

对于表面结构，有以下几个需要注意的:

  • xy 方向上表面的大小（15 Å或者100原子左右）;
    A 这个影响表面吸附物种的覆盖度;
    B 影响体系的尺寸大小和计算时间;
    C 不同的大小需要选取对应的K点；回想一下我们前面提到的经验规则。
  • 不同的晶面，(111)， (100)，(110);
    A 这取决于你研究的方向;低指数面、解离面。
    B 不同晶面的表面能;
    C 不同晶面的表面结构，反应活性等。
  • 表面结构的层数（3+3以上，6Å+6Å）
    A 层数多了，原子变多，体系在 z 方向的尺寸增加，也会影响计算速度；
    B 计算中需要弛豫的层数；
    C 不同层数对你要计算的性质会产生影响，比如表面能；
    D 不同晶面需要的层数也不尽相同，一般开放的表面需要更多的层数；
    E 根据自己感兴趣的性质，选择合适的层数，也就是需要你去测试一番。
  • Slab模型有两种，一种是上下表面对称的，一种是非对称的。对称性结构往往需要很多层，体系较大。 非对称的结构体系较小，但存在偶极矩的影响，要注意加LDIPOL 和IDIPOL这两个参数来消除。1、2、3分表代表在x、y和z 方向上进行校正，4代表在所有方向上。真空层15-20埃。

3. 工具：MS、ASE、P4vasp、vaspkit。

MS操作。

一般用conventional cell。
MS操作。Windows
VESTA将CONTCAR转成cif,导入MS。
build-surface-cleave surface 切面
Surface mesh-surface vectors 可调整vector。

Surface box-top/thickness，调整终端和厚度。
Build-crystals-build vaccum slab-图2-build。

Lattice parameters-advanced-re-orient to startard把真空层沿Z轴 图3。

转换成POSCAR，可优化后进一步扩胞。

ASE操作。python

块体优化后
module load python/3.7.3
python

from ase.io import read
from ase.io import write
from ase.build import surface
a=read(”CONTCAR”,format=’vasp’)
b=surface(a,(1,1,1),6,vacuum=7.5) 6层，真空层15Å
write(”POSCAR”,b,format=’vasp’)
Ctrl+D （退出 python）

大师兄ASE脚本：切金属

python3 cssm.py

vaspkit

块体优化后复制成POSCAR
vaspkit –task 803
生成SLAB***.vasp

4 表面弛豫计算

单点能

转换成POSCAR单点计算。或表面优化时第一个离子步的能量，但要保证这一步收敛。
KPOINTS设置，20-30。
INCAR设置，NSW=0或默认0.偶极矩矫正 LDIPOL = .TRUE. ; IDIPOL = 3。
POTCAR准备，提交任务脚本。完成后 grep without OUTCAR 获取单点能。

表面优化

1）POSCAR第7行后面加一行加S， Sed –i ‘7a S’ POSCAR
原子坐标加F/T。 Sed –i ’10,11s/$/ F F F/g’ POSCAR
Sed –i ’12,13,14s/$/ T T T/g’ POSCAR
vaspkit –task 402/403 或者ctrl+v,选中列-大写A-输入 T T T - ESC.
2）POTCAR和KPOINTS不变。
INCAR设置：IBRION=1、2，EDIFF=1E-6，EDIFFG=-0.03（比块体小一点），ISIF=2，
提交任务。

5 表面能计算

表面能是创造物质表面时，破坏分子间化学键所需消耗的能量。在固体物理理论中，表面原子比物质内部的原子具有更多的能量，因此，根据能量最低原理，原子会自发的趋于物质内部而不是表面。表面能的另一种定义是，材料表面相对于材料内部所多出的能量。把一个固体材料分解成小块需要破坏它内部的化学键，所以需要消耗能量。如果这个分解的过程是可逆的，那么把材料分解成小块所需要的能量就和小块材料表面所增加的能量相等。但事实上，只有在真空中刚刚形成的表面才符合上述能量守恒。因为新形成的表面是非常不稳定的，它们通过表面原子重组和相互间的反应，或者对周围其他分子或原子的吸附，从而使表面能量降低。
γ=($E_{surf}$ - N * $E_{bulkatom}$)/(2A)

1ev/$Å^2$ =16.02J/$m^2$
$E_{surf}$是未弛豫的表面能量,A是表面积a*b。

吸附能扫描

算吸附能时一般只需要算分子在高对称点的吸附能，比如算Li在铜表面的吸附能，我们只需算Li在Cu(111）面的Top，HCP，FCC和Bridge位的吸附能即可。但是这些离散的点没法帮你构建一个“势能面”，也就是无法得到吸附能在表面究竟是怎么分布的. 下面的文献里DOI：10.1063/1.4901055有一张图非常生动，讲的分别是Li在Li(001)以及Mg在Mg(0001)表面的吸附势能分布，我们不仅可以知道哪些位置有利于吸附，还可以根据吸附能分布，描绘出分子在表面的扩散路径，这有助于我们使用NEB方法算扩散能垒。

图片

其实思路很简单，就是将表面网格化，算分子在格点的吸附能，再画出等高线图。但是大部分吸附位点其实是不稳定的，所以我们采取固定吸附分子x,y方向，在z方向弛豫以达到吸附平衡的策略。

其实手动撒点，再采集数据也是可行的，但是会比较麻烦，因此作者根据实际需要开发了一款脚本scan_adsorption_energy用于自动完成这个过程。脚本使用Python编写，需要numpy和matplotlib第三方库。 我们首先算好一个吸附例子得到CONTCAR,这个可以让我们得到吸附分子的元素信息和理想的吸附高度。如下图，吸附的Li是第97号，也是最后一个原子，接着就可以变换这个结构得到不同的吸附结构。

图片

首先我们需要在目录下准备以下几个文件CONTCAR,INCAR,KPOINTS,POTCAR以及排队脚本文件(我的是vasp.pbs)，然后使用脚本产生不同的结构：

tamas@tamas-PC:~/Desktop/scan$ python scan_adsorption_energy.py -g
Input which atom you want to move to scan adsorption energy–>-1
脚本提示我们移动第几个原子产生不同的吸附结构，我们可以输入97或者-1，移动最后一个Li。这个Cu(111)面是(331)的超胞，因此我们最需要取左下角的1/9的小区域撒点即可，然后程序会将小区域内的数据点进行平移得到整个表面的数据。如何根据不同的表面调整这个参数以及撒点的密度呢？可以使用-u参数进行自定义设置。

Input how many supercells in x axis–>3
Input how many supercells in y axis–>3
Input interporate how many points for each unit cell in x axis,even number are suggested–>2
Input interporate how many points for each unit cell in y axis,even number are suggested–>2
command to submit jobs,e.g. qsub vasp.pbs–>qsub vasp.pbs
Input which atom you want to move to scan adsorption energy–>-1
手动选择x方向的超胞数3，y方向的超胞数3,在小区域内的x方向的撒点数，y方向的撒点数（为了包含中心点，建议选择偶数）和提交排队脚本的命令（我的是qsub vasp.pbs）。脚本会自动产生4个结构（2*2），并会将每个计算所需的文件都复制到相应的文件夹内，并提交排队脚本进行计算。同时脚本会写出一个grid.info里面包含了超胞数，选择移动的原子以及晶格常数信息。 在所有的计算完成后，同样使用该脚本提取数据。

python scan_adsorption_energy.py -e
脚本会自动从每个结构文件夹内读取能量数据和吸附原子的坐标。生成两个dat文件new_direct.dat和new_dat-3D.dat。前者是原子的分数坐标和总能量，后者有4列，分别为原子的笛卡尔x分量坐标，y分量坐标，z分量坐标以及总能量。我们需要使用第1,2,4列绘制等高线图。如果需要绘制吸附能，需要自己手动将第四列的数据减去吸附原子和表面的能量。 脚本会尝试使用matplotlib绘制等高线图，如下图，你也可以使用Origin Pro绘图。

图片

脚本可以在我的Github仓库下载(https://github.com/tamaswells/VASP_script/blob/master/scan_adsorption_energy.py).

固定基矢优化结构

在第一性原理计算时，二维材料的建模会采用在某一方向（一般是z轴方向）添加15-20 Å的真空层的方式来抵消层与层之间的相互作用。但是，当我们优化晶格常数的时候会遇到一个问题：那就是使用ISIF=3来优化平面内晶格常数（一般是a和b）时，c轴方向的晶格长度总是随着优化变的越来越小。

  1. 采用ISIF=4的固定晶胞体积的方法来计算。这样采用实验参数建的晶胞，由于a和b方向的变化，c轴也会变化。解决了越是优化c轴越小的问题，但是计算量还是很大。

  2. 编译vasp，就自己编译了固定z轴的脚本，专门计算二维材料。这种方法优点是直接换脚本就能解决上面优化时c轴缩小的问题，但是我想固定两个轴就需要再编译一份脚本。比较麻烦。

修改constr_cell_relax.F实现固定基矢优化结构的原理：

FCELL contains the forces on the basis vectors. These forces are used to modify the basis vectors according to the following equations:

FCELL 是个3 * 3的矩阵，包含了晶格矢量的受力，当ISIF = 3时，晶格矢量会按照如下的方法更新：

A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
DO J=1,3
DO I=1,3
DO K=1,3
A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
ENDDO
ENDDO
ENDDO
当通过OPTCELL文件把晶格矢量在某个方向的受力改成0以后，

IF (ICELL(I,J)==0) FCELL(I,J)=0.0

在结构更新的时候，该方向就不会改变。完成编译后，在 INCAR 文件中设置 ISIF=3，便可实现固定 Z 轴的结构优化。
对于其他方向晶格基矢的修改同理：对于a方向基矢，将 FCELL(3,I) 和 FCELL(I,3) 分别改为 FCELL(1,I) 和 FCELL(I,1) ；对于b方向基矢，则分别改为 FCELL(2,I) 和 FCELL(I,2) ；固定两个基矢则应该同时将两个方向对应的矩阵元设置为0。或参考 VASP并行可执行软件包，可对晶胞参数进行部分优化 可实现添加参数文件实现不同基矢方向固定。（未做测试）

  1. 肖承诚博士的补丁，按照肖博士的方法安装以后，在INCAR中添加如下参数即可。

IOPTCELL = xx xy xz yx yy yz zx zy zz

上面的九个参数对应晶格常数矩阵。

1 xx xy xz

2 yx yy yz

3 zx zy zz

1代表驰豫，0代表固定

此时固定c轴的设置就是如下：

IOPTCELL = 1 1 0 1 1 0 0 0 0

下载地址：https://github.com/Chengcheng-Xiao/VASP_OPT_AXIS

滑移能和剥离能

材料的塑性行为主要是由于材料的结构和化学键决定的。材料发生塑性形
变在原子尺度上可以分解成两个过程。

一个是由于原子、位错或者晶界沿着某
一特定方向自由运动；第二个是在这种运动过程中存在束缚力保证在其它方向
不会产生裂纹，导致材料的断裂。因此任何一个塑性材料必须满足两个要求，
一是对某些原子层存在势垒相对低的滑动方式，以确保原子、位错或者晶界可
以沿着这些方向滑动。另一个是在这些滑动过程中，滑移平面两侧的原子依然
存在某种相对强的作用力，以确保滑移过程中材料的完整性。如金属键，其自
由电子容易发生变形和转移，在滑移过程中金属键不断的断裂和产生，从而导致了金属的滑移能量很小。而金属键的键强度很大，保证的金属滑移过程的完
整性。而离子键和共价键由于配位数和方向性等要求，在滑移时必须先发生键
的断裂，从而阻碍滑移的进行。由范德华力连接的层状材料，虽然层与层之间
容易发生滑移，但其层间作用力较弱，材料在滑移过程中容易发生破坏。

为了研究 α-Ag2S 的塑性行为，我们通过第一性原理对 α-Ag2S 的滑移过程进
行了计算。其中，我们计算的滑移面为（100）方向，滑移方向为[001]方向。其
滑移过程如图 6.10 所示，α-Ag2S 的滑移面沿着能量最低点进行滑移，滑移面不
仅发生相对移动，其面间距也会发生。因此，整个滑移路径成波浪形。 我们将滑移一个周期（滑移一个晶格常数）均分 12 步。每一步我们计算系
统总能量与滑移面之间距离的关系，对于塑性材料而
言，它需要同时具有足够小的滑移能来保证材料的滑移可以进行；以及大的解离
能以至于材料不会发生断裂。