VASP通过CI-NEB计算气体在孔内的扩散

突然发现我的博客里面竟然没有我用的最多的vasp的教程，这怎么可以，正好最近在开展一个关于气体扩散的课题，我就想着用CI-NEB来计算一下气体的扩散能垒，用以验证实验结果的推论。

前言

这次课题是一个炭材料脱氢的工作，大概就是需要我们计算一下氢气在特定的炭材料孔洞中的扩散能垒，具体的工作中肯定是要计算好多个模型来进行对比的，我这里作为一个教程，就只展示一个模型的计算了。

关于CI-NEB的计算其实我之前做了这么久的vasp相关工作，但基本都没有涉及过，这次也算是首次尝试了，有不对的地方以后就再回头修改了，主要还是参考大师兄的那一套，不过很幸运的是我们课题组的vasp是编译了vtst的，所以可以直接省略不少前期工作。以下是大师兄的网址：

Ex70 过渡态的计算（一） | Learn VASP The Hard Way (bigbrosci.com)

大师兄的博客把整个流程分了好多章节，而我一个拖延症患者，我还是觉得一次性写完对我以及对需要这样课程的人来说都是比较好的事儿。

前期准备

CI-NEB的具体原理我就不多做赘述了，引用大师兄博客里面的一段话吧，然后需要注意的是你需要vtst这个脚本包，它的获取还是很方便的，可以自己Google或者百度。

算过渡态，肯定指的是某一个反应或者过程的过渡态。过渡态的两边也分别对应着反应物和反应产物。所以，算一个反应的过渡态，不管计算过程如何，我们最终都会有三个结构：初始和末态结构，以及过渡态。

这个过程中，我们需要先获得初态和末态，这个就是根据实验情况来进行建模了，以下是我建立的模型：

对这两个模型进行优化，但是其实作为周期性结构，这两个模型中的氢气分子实际上距离是非常近的，也就是说，如果你用这两个模型来做出初态和末态进行插点的话，会导致它的中间态位移非常小，且不会经过中间官能团的区域，而是从初态出发，往左移动，越过盒子的左边界，从末态的右边界出现，再逐渐到末态的位置。所以这样的模型是很不合理的，但是我又不想只用一半的模型，毕竟这里三个含氧官能团呢，所以我们可以在这两个状态之间加上一个mid-state，这样扩散路径就可以被指定为ini-mid，mid-fin，两步。以下是我的mid-state的模型：

所以在最初，我们需要对这三个模型进行建模，然后对它们进行优化，这里KPOINTS我用的是221，以下是我结构优化所用的INCAR文件：

SYSTEM = H2

# === Basic Parameters ===
LWAVE        = .False
LCHARG       = .False
PREC         = Accurate
ISTART       = 0

# === Electronic Relaxation ===
EDIFF = 1E-5
ENCUT = 400
LREAL = Auto

# === Ionic Relaxation ===
EDIFFG = -0.02
IBRION = 1
ISIF = 2
NSW = 500

# === DOS related parameters ===
ISMEAR = 0
SIGMA = 0.02

# === DFT-D3 ===
IVDW = 11

NCORE = 32

优化结束之后我们得到这三个模型的CONTCAR，那么我们就可以继续下一步了。

ini-mid与mid-fin的工作内容其实是一致的，其实如果自己有时间倒腾的话，甚至可以放弃mid-state，插点的过程自己用手去调整。

将ini和mid的CONTCAR复制到同一个目录下，注意重命名一下，我们能可以将它们命名为ini-POSCAR和mid-POSCAR，然后我们这里需要配置一下vtst，其实就是添加一个环境变量的工作，我这里下载的是 vtstscripts-1022版本，然后把它保存到了我自己存放脚本的目录下，然后进到vtstscripts-1022，先给我们需要用到的插点的脚本，也就是nebmake.pl施加可运行权限，在命令行输入：

chmod +x nebmake.pl

这里可以先用pwd指令查看一下当前目录的路径，便于后面修改环境变量，然后输入：

# 查看当前路径
pwd
# 系统返回的路径
/home/jason/software/script_file/vtstscripts-1022
# 配置环境变量
vim ~/.bashrc

在最后一行添加：

export PATH="/home/jason/software/script_file/vtstscripts-1022:$PATH"

可以看到中间的路径就是我们刚刚pwd得到的路径，然后再在命令行中输入：

source ~/.bashrc

这个时候nebmake.pl就可以在任意目录下正常使用了

插点并提交任务

在上一部分我们的工作目录下已经有了ini-POSCAR和mid-POSCAR，并且配置好了nebmake.pl，那么我们即可以开始插点了，这里插点的数量不是越多越好，适宜即可，所以我们选择插入四个点，这样算上mid-fin的过程，我们整个过程插入了八个点，再加上三个基态的点，总共就会有11个点。

接下来说说用nebmake.pl插入点，这个脚本的使用方法如下：

# 第一个为初态的结构文件，第二个为末态的结构文件，第三个为插入的点数量
nebmake.pl IS FS N

所以我们要插入四个点，那么就可以用以下的指令：

nebmake.pl ini-POSCAR mid-POSCAR 4

使用了这个指令后目录下会多出来00，01，02，03，04，05六个文件夹，每个文件夹都分别保存了一个POSCAR，其中00和05目录下的POSCAR就是我们刚刚输入的ini-POSCAR和mid-POSCAR，而另外四个文件夹就是nebmake.pl帮我们生成的插点POSCAR，注意这个插点不一定是可靠的，我们一定要自己用可视化软件打开这些插点的模型查看一下。

这就是我生成的四个点，但是其中第一个和第二个模型由于生成的点距离下层的羟基有点近了，所以我在vesta中对它们的坐标进行了一些小调整，那就是把它们的z轴坐标适当增大了一些，让氢气不至于离羟基太久而导致成键（成键会导致能垒变得很大，计算的过渡态也就没有啥意义了，我们这只是氢气的扩散，不涉及成键过程）

检查没有问题之后我们就可以开始准备提交任务了。这个时候我们的工作目录下理论上来说是有六个文件夹以及两个POSCAR文件，即ini-POSCAR和mid-POSCAR文件，这俩文件在插点完成之后其实就没啥用了，可以删去也可以放在这里不动它，我们还需要在工作目录下准备计算的其它文件，也就是另外的三件套，INCAR，POTCAR以及KPOINTS，K点我们依然设置为221，POTCAR可以沿用之前结构优化的POTCAR，但是INCAR就需要修改一下了，一下是我用的INCAR：

SYSTEM = H2

# === Basic Parameters ===
LWAVE        = .False
LCHARG       = .False
PREC         = Accurate
ISTART       = 0

# === Electronic Relaxation ===
ENCUT = 400
LREAL = Auto

# === DOS related parameters ===
ISMEAR = 0
SIGMA = 0.02

# === DFT-D3 ===
IVDW = 11

# === CI-NEB ===
ISIF = 2
LCLIMB = .TRUE.  # 用CI-NEB方法算TS
IMAGES = 4       # 反应路径插了4个点
IBRION = 1       # 初始结构（路径）感觉不太好，可以用 IBRION = 3
POTIM = 0.2      # 一般0.2可以满足大部分的情况
SPRING = -5      # 约束不同IMAGES之间的力
EDIFF = 1E-4     # 电子步的收敛标准
EDIFFG = -0.05   # 离子步的收敛标准，在粗算中就是个摆设
NSW = 60         # 粗算中50-100均可

NCORE = 32

第一步的粗算，主要还是看看我们插入的点是否合理，结构优化之后是否会散掉或者变形。

然后就可以提交任务开始计算了。

数据的后处理

经过一段时间的计算，计算完成，就可以对计算结果进行处理了，这里首先最重要的就是查看能量的变化，查看能量的变化理论上来说vtst自带的脚本就可以实现，但是不知道为啥我这个不好用，所以我还是自己参考大师兄的教程写的一个小脚本：

#!/bin/bash
for i in *; do
if [ -e $i/OUTCAR ]; then
echo -e  $i "\t" $(grep '  without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
fi
done

将这些内容保存为energy.sh，然后直接在工作目录下输入：

sh energy.sh

运行脚本即可，就会自动将所有插入的点的能量输出，然后就可以输出数据进行处理了。