0%

算过渡态,光知道检查虚频是不对的。如果算出来虚频多,很大程度是因为你的NEB初始结构有问题,也就是说你的过渡态路径背后的物理化学意义不是那么地理想。而对于NEB的初始结构,大部分人都是通过VTST的nebmake.pl脚本来实现的。

1
nebmake.pl IS FS 8 

分析上面这个命令:有4个关键的信息:

  • 脚本:nebmake.pl
  • 初始结构: IS
  • 末态结构: FS
  • 插的点数:8

本节默认你优化好了IS,FS,并且插8个点,分析一下脚本nebmake.pl的2个坑。其实,在Density Functional Theory: A Practical Introduction 这本书的第6章,其中的一个坑已经讲过了。强烈建议新手,老手认真看看这一章。

想知道这两个坑,首先我们先分析下nebmake.pl的工作原理。简举个化的一维例子:一个线段两端的坐标是x1和x2,我们把这个线段分成n份,并获得每一段的起始点的坐标。想必大家都知道怎么弄。

x_i = x0 + i *(x2-x1)/(n+1)

如果扩展到三维的xyz坐标,分别对y、z进行同样的处理。我们就得到了初始结构和末态结构之间这些的IMAGES的坐标。具体见:第六章的148页。

知道了原理,现在就可以分析其中的两个坑了。

坑1

前面讲的原子在表面扩散的例子,如果原子从fcc扩散到hcp。fcc和hcp的结构中,原子在z方向的坐标基本是相同的。如果运行nebmake.pl 这个命令,就会出现这样的一种情况,所有的IMAGE结构中,z的坐标几乎不变,也就是原子在表面横着走。而实际情况呢?原子在扩散的过程中,z方向的坐标也发生相应的变化。就如同你从山的一边爬到另一边,虽然海拔基本没变,但爬山这个过程,有升,下山过程,有降。所以,这种情况下,直接用nebmake.pl插的点在z方向上的物理意义并不准确。稍微扩展下,在x或者y方向上也可能会发生类似的情况,如果IS和FS在某一维度的变化很小时,一定要注意这一维度上的物理意义是否可以被表现出来。

既然知道了这个坑,我们该怎么填呢?

方法1:

管它是不是坑,或者不知道这个坑,直接开车冲过去。这是很多新手常见的做法,虽然有时候可以开过去,但不推荐,毕竟也会溅一车泥,搞不好坑大了还会掉进去。

方法2:

学习书中的例子,在计算之前先手动修改下结构。比如插了8个点,把第2和第7个IMAGES中的坐标向上移动0.1$\AA$。第3和6个向上移动0.15$\AA$,第4和5移动0.2$\AA$。这样做的好处就是,初始结构在z方向上具有更好的物理意义,使得计算收敛更快。

上图是书中的一个例子,最上面的曲线是没有调节z坐标时,初始结构所对应的能量;中间的为调节初始结构的z坐标后,初始结构的能量;最下面的是NEB计算完成之后,各个IMAGE的能量。从能量上也可以看出来,如果z坐标我们不修改的话,体系能量与稳定的相差甚远,间接告诉我们可能需要更多的优化步数来收敛。这个强烈建议大家自己亲手算一算。体系简单,不耗费那么多机时,有助于加深对NEB的认识和学习。

前面几节,我们从gamma点到3x3x1计算过渡态。首先复习一下计算的流程:

1) 我们计算的是金属表面上H原子的扩散;

2) 使用gamma点计算的时候,slab固定住了;

3)gamma点计算结束后,检查结果:

  • 看到有一个漂亮的NEB图;
  • 能量变化也是稳稳妥妥滴;
  • 结构查看也没啥异常情况;

4) 在确认第3)步检查OK之后,将gamma点的计算备份;增大K点至3x3x1继续算。

5)计算结束后,重复第3)步的检查,确认没啥问题。

继续算

前面这么做很啰嗦(大师兄本人也很啰嗦),目的只有一个:用最少的机时获取最好的NEB初始结构。当我们完成这一步之后,就可以再继续下面的操作:

  • 备份3x3x1的计算;

  • 增大K点至5x5x1;

  • 放开表面的原子;(POSCARtoolkit.py
  • 继续算,结果如下:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ ls
00 02 04 06 08 INCAR POTCAR job_sub movie neb.dat slurm-1133307.out spline.dat vasprun.xml 01 03 05 07 09 KPOINTS exts.dat mep.eps movie.vasp nebef.dat slurm-1133315.out vaspgr
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ cat KPOINTS
K-POINTS
0
Gamma
5 5 1
0 0 0
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OUTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OUTCAR
User time (sec): 3864.268
System time (sec): 15.316
Elapsed time (sec): 3905.371

Maximum memory used (kb): 466776.
Average memory used (kb): 0.

Minor page faults: 816554
Major page faults: 29
Voluntary context switches: 29864
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OSZICAR -n 1
10 F= -.32111627E+03 E0= -.32109808E+03 d E =-.130617E-03
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ ta.sh
00 -321.15187450
01 -321.13494327
02 -321.09807672
03 -321.06314408
04 -321.03855868
05 -321.02460382
06 -321.01884134
07 -321.02967204
08 -321.06853536
09 -321.09850728

从上面可以看到,顶点在06的位置,差不多就是过渡态了。

验证过渡态

那么我们算出来的过渡态到底对不对呢?

上面,我们通过能量分析,06这个Image就是过渡态了;但06就真的是过渡态吗?下面我们需要做2件事情:

1) 查看结构:

为了区分明显,上图中桥式位置两端的Ru原子,用暗红色标记出来。可以看出,06结构中,H原子在桥式的位置上;是我们想要的过渡态。结构这一关也过了。

2)频率分析:对于一个基元反应的过渡态来说,会有一个对应的虚频。因此,我们将06的CONTCAR单独取出来,做一个频率分析:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ mkdir freq && cp 06/CONTCAR freq/POSCAR && cd freq &&sed -i '10,27s/T/F/g' POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ cp ../INCAR .
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ vi INCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ kpoints.sh 1 1 1
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ pospot.sh
Generating NEW POTCAR...
************************
Done
************************
NEW POTCAR containes....
Ru H
************************
Elements in POSCAR
************************
Ru H
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ qsub
Submitted batch job 1179250
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ ls
CHG CONTCAR DYNMAT IBZKPT KPOINTS OUTCAR POSCAR REPORT XDATCAR slurm-1179250.out CHGCAR DOSCAR EIGENVAL INCAR OSZICAR PCDAT POTCAR WAVECAR job_sub vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ grep cm-1 OUTCAR
1 f = 36.804765 THz 231.251161 2PiTHz 1227.674788 cm-1 152.212329 meV
2 f = 30.534675 THz 191.855022 2PiTHz 1018.527097 cm-1 126.281311 meV
3 f/i= 9.822527 THz 61.716759 2PiTHz 327.644231 cm-1 40.622722 meV
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$
  • 创建freq文件夹,将06的CONTCAR复制到freq中;
  • 将表面的Ru原子固定住;
  • 复制NEB计算的INCAR过来,修改频率计算相应的参数;
  • 使用gamma点算频率;
  • 生成POTCAR;
  • 提交任务。
  • 查看结果,有一个虚频。计算完事!

注意:

  • 不会算的看前面频率计算的内容;

  • 自己测试下不同K点:2x2x1;3x3x1; 5x5x1计算出来的零点能有什么区别?频率大小有什么区别。

小结:

拖拖拉拉,磨磨唧唧,终于把H原子在表面过渡态计算的NEB讲的差不多了。如果你一路沿着教程走过来,最终看到一个虚频对的结果时,会发自内心由衷地一笑:啊!啊!啊!过渡态原来这么简单。为避免这种情况的发生,我先给你脑勺拍一板凳,后面事还多着呢。还是那句话:做任何一个事情,从你认为它简单的那一刻起,你就输了。算过渡态,

1)要把初始结构,末态结构优化好;

2)合理采用粗糙的模型,来检查自己设想的反应路径;

  • 粗糙的模型:一方面指的是体系的大小,一方面是计算参数的设置;
  • 切忌直接上来硬算,要不然把服务器累个半死,却得不到多少好的结果。

3) 从能量,结构,虚频等多个角度去分析你的结构。

本节计算的文件已经打包上传到百度网盘,由于版权,压缩包里面的POTCAR就不放了,大家自己生成一下。下载链接:链接:https://pan.baidu.com/s/1sfOLxTB5Rdr5Il7vNCOkyA 提取码:m591

什么是ASE?

ASE 是Atomic Simulation Environment的简称。详细的介绍大家可以浏览官网:https://wiki.fysik.dtu.dk/ase/about.html。ASE是一个基于python的程序库,可以做的事情很多,不同程序任务的设置,结构的搭建,结构分析都可以做。

为什么用ASE?

  • ASE里面包含了很多的功能,其中个人认为最为方便地就是计算模型不同文件类型的转换。比如将cif文件转化为POSCAR,将POSCAR转化为xyz文件,等等;

  • 可以在学习ASE的时候,顺便学习python的一些知识;

  • 使用ASE分析计算结果;

  • 使用ASE学习做计算;这里分享江山推荐的一本书:https://github.com/jkitchin/dft-book

怎么安装ASE?

如果你用的是国防科大吕梁超算中心,管理员已经安装好了ASE,那么需要你做的只有一步,添加ASE的路径到.bashrc文件中就可以直接用了。

1
2
3
4
5
6
iciq-lq@ln3:/THFS/opt/python3.6/bin$ ls ase*
ase ase-build ase-db ase-gui ase-info ase-run
iciq-lq@ln3:/THFS/opt/python3.6/bin$ pwd
/THFS/opt/python3.6/bin
iciq-lq@ln3:/THFS/opt/python3.6/bin$ tail -n 1 ~/.bashrc
export PATH=$PATH:/THFS/opt/python3.6/bin

如果你想在自己电脑上安装ASE。安装过程很简单,详细的信息请参考:https://wiki.fysik.dtu.dk/ase/install.html。Mac太贵,没尝试过,Windows系统本人不太熟悉,下面附上在Windows10中Ubuntu寄生系统的安装过程(跟正儿八经Ubuntu系统的安装流程一样); 总共就三步。

  • 安装pip:如果已经安装了pip,则就剩2步了。

    1
    2
    sudo apt install python-pip
    sudo apt install python3-pip

    一个是基于python2,一个是python3的。

  • 安装ASE:

    1
    pip install --upgrade --user ase

  • 添加ASE的路径到.bashrc文件。

1
export PATH=$PATH:/home/bigbro/.local/bin

在Terminal里面输入ase,然后摁下tab键:如果跟下面的输出一样,说明八九不离十了。

1
2
3
4
bigbro@DESKTOP-S6R3TUP:~$ ase
ase ase-build ase-db ase-gui ase-info ase-run ase_figures.sh aselite.py aselite.pyc
bigbro@DESKTOP-S6R3TUP:~$ ase --version
ase-3.17.0

前面一节我们通过粗糙的方法演示了一遍怎么在超算上跑过渡态的NEB计算。一般来说,通过这一步我们可以大体上确定这个基元反应中的路径是不是靠谱的了。不论你得到一个粗糙的NEB能量曲线:

亦或者一个比较精细的:

这都表明,你快要得到这个过渡态结构了。在进行后面操作讲解之前,我们先做2个事情:

1) 由于本人也是老司机了,算过渡态的时候遇到的那些零零碎碎的问题一时也很难记起来,如果你正在算过渡态,并且遇到了问题和错误,可以将计算打包,以及遇到的问题发给我。如果本人有幸可以帮你解决问题,那么,作为交换条件,这个例子就写在后续的教程里面。有下面几点需要注意的:

  • i)请使用邮件发送例子(lqcata@gmail.com),尽可能保证计算的文件完整,千万不要删掉某些文件后再打包发给我(CHG,CHGCAR,WAVECAR除外)。如果例子太大,可以上传到百度网盘,然后将链接和你遇到的问题发给我。

  • ii)本人最近很忙,可能回复有些慢,请见谅。如果你看到我在QQ群里面瞎BB的话,可以提醒下我去查看邮件。

  • iii)只想解决自己问题,不愿意分享自己例子的,请不要发邮件给我。

2)第二件事就是分享一个自己平时常用的脚本,这个脚本主要是将前面的计算步骤保存起来,然后将CONTCAR复制成POSCAR用以下一步的计算。脚本的名字为:save_calculations.sh.

1
2
3
4
5
6
7
8
9
10
#!/usr/bin/env bash 

mv POSCAR POSCAR-$1
mv OUTCAR OUTCAR-$1
mv OSZICAR OSZICAR-$1
mv vasprun.xml vasprun.xml-$1
mv EIGENVAL EIGENVAL-$1
mv IBZKPT IBZKPT-$1
cp CONTCAR CONTCAR-$1
mv CONTCAR POSCAR

脚本演示:

1
2
3
4
5
6
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ ls
CONTCAR DOSCAR DYNMAT IBZKPT INCAR KPOINTS OSZICAR OUTCAR p4vasp.log POSCAR POTCAR sub12 sub24 vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ save_calculations.sh gamma
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ ls
DOSCAR DYNMAT IBZKPT INCAR KPOINTS OSZICAR-gamma OUTCAR-gamma p4vasp.log POSCAR POSCAR-gamma POTCAR sub12 sub24 vasprun.xml-gamma

脚本说明:

  • 脚本内容很粗暴直接,大家嫌不好看也可以结合下for循环改进地更加简洁些,本人也懒得改,能实现目的就OK了;
  • 效果很简单,就是把之前计算的一些结果重新命名一下,大家可以根据自己的爱好或者习惯,将gamma换成其他的。

提交任务

我们把参数设置地稍微精细一下,然后再跑一遍NEB。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls *
INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml

00:
POSCAR

01:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT WAVECAR XDATCAR
.
.
.
08:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

09:
POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ for i in 0{1..8}; do cd $i ; save_calculations.sh gamma ; cd - ; done
mv: cannot stat `vasprun.xml': No such file or directory
/THFS/home/iciq-lq/LVASPTHW/ex78
.
.
.
mv: cannot stat `vasprun.xml': No such file or directory
/THFS/home/iciq-lq/LVASPTHW/ex78
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls *
INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml

00:
POSCAR

01:
CHG CHGCAR CONTCAR-gamma DOSCAR EIGENVAL-gamma IBZKPT-gamma OSZICAR-gamma OUTCAR-gamma PCDAT POSCAR POSCAR-gamma REPORT WAVECAR XDATCAR
.
.
.
08:
CHG CHGCAR CONTCAR-gamma DOSCAR EIGENVAL-gamma IBZKPT-gamma OSZICAR-gamma OUTCAR-gamma PCDAT POSCAR POSCAR-gamma REPORT stdout WAVECAR XDATCAR

09:
POSCAR

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ kpoints.sh 3 3 1
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ yhbatch -p gsc -N 2 -J test job_sub
Submitted batch job 1133307
  • 这里我们把01到08中的一些输出文件都保存了一下,然后将CONTCAR复制成POSCAR用来进行下一步的计算;
  • 运行的时候,会报错mv: cannot stat ‘vasprun.xml’: No such file or directory,不用管就行,因为images的文件夹中没有vasprun.xml`文件,嫌报错烦的可以自己加个if语句修改下脚本。
  • 我们将K点增加到331,继续计算。

查看计算结果

查看结果的时候,有以下三个主要的方面:

1) 看能量

  • 可以使用VTST的nebresults.pl 的小脚本,也可以使用自己写的小脚本(ta.sh)。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 01 02 03 04 05 06 07 08 09 exts.dat INCAR job_sub KPOINTS mep.eps movie movie.vasp neb.dat nebef.dat POTCAR slurm-1133307.out spline.dat vaspgr vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.023782 -321.151800 0.000000
1 0.007650 -321.196900 -0.045100
2 0.009671 -321.160900 -0.009100
3 0.010446 -321.128300 0.023500
4 0.014498 -321.106600 0.045200
5 0.019054 -321.096100 0.055700
6 0.019424 -321.094300 0.057500
7 0.013325 -321.111700 0.040100
8 0.005225 -321.156900 -0.005100
9 0.011589 -321.098500 0.053300

Extremum 1 found at image 0.910082 with energy: -0.046551
Extremum 2 found at image 5.724723 with energy: 0.057913
Extremum 3 found at image 8.102362 with energy: -0.007602
Extremum 4 found at image 8.999187 with energy: 0.053367

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ta.sh
00 -321.15187450
01 -321.19698136
02 -321.16098546
03 -321.12832705
04 -321.10664059
05 -321.09613563
06 -321.09434359
07 -321.11177677
08 -321.15691232
09 -321.09850728

  • 看到能量的时候,前面我们粗算的时候,已经知道能量随着路径的变化是怎么样的,这里我们就脑子自动画图,从01到08,看着能量慢慢上升,在05,06左右的时候达到最高点,然后继续下降。

  • 注意:我们这是接着前面一节进行的计算, 00 和 09中的能量都是551K点下的能量,因此得到的能量随路径的变化是下面这样子的。如果看到这样的图片,不要慌张,跳过第一个和最后一个点,直接看中间的部分即可。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep irre */OUTCAR
    00/OUTCAR: Found 9 irreducible k-points:
    01/OUTCAR: Found 5 irreducible k-points:
    02/OUTCAR: Found 5 irreducible k-points:
    03/OUTCAR: Found 4 irreducible k-points:
    04/OUTCAR: Found 4 irreducible k-points:
    05/OUTCAR: Found 4 irreducible k-points:
    06/OUTCAR: Found 4 irreducible k-points:
    07/OUTCAR: Found 4 irreducible k-points:
    08/OUTCAR: Found 5 irreducible k-points:
    09/OUTCAR: Found 13 irreducible k-points:

2) 看计算有没有正常结束。

主要是查看:i) NEB跑了多少步,ii) 自己设置了多少步,iii) 通过OUTCAR中结构收敛的特定关键词。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 02 04 06 08 exts.dat job_sub mep.eps movie neb.dat POTCAR spline.dat vasprun.xml
01 03 05 07 09 INCAR KPOINTS mep.png movie.vasp nebef.dat slurm-1133307.out vaspgr
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ tail 03/OSZICAR -n 5
RMM: 2 -0.318010871323E+03 -0.48063E-06 -0.15546E-06 134 0.370E-03 0.397E-03
RMM: 3 -0.318010871353E+03 -0.29988E-07 -0.12836E-08 21 0.273E-03 0.402E-03
RMM: 4 -0.318010871104E+03 0.24867E-06 0.00000E+00 17 0.260E-03 0.157E-03
RMM: 5 -0.318010871201E+03 -0.96799E-07 -0.11232E-08 18 0.266E-03
14 F= -.32114179E+03 E0= -.32112833E+03 d E =-.102536E-03
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep NSW INCAR
NSW = 200
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep 'reached re' 03/OUTCAR
reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$

当然你可以写一个小脚本如下来避免每次都重读在终端里面打出来。

1
2
3
4
5
#!/usr/bin/env bash 

tail 03/OSZICAR -n 5
grep NSW INCAR
grep 'reached re' 03/OUTCAR

也可以设置下alias来判断:

1
alias check='tail 03/OSZICAR -n 5 && grep NSW INCAR && grep "reached re" 03/OUTCAR '

有能力的可以写一些更高级点的,加些if,for语句智能判断NEB有没有算完,有没有收敛。

3) 查看NEB中各个IMAGE中的结构有没有跑乱,散架,有没有物理化学意义。前面我们讲过了几种办法,本节我们简单复习一下,就不再啰嗦了。

i) 通过vaspkit 结合 vmd实现看动画的效果(Windows用户)

ii)使用ASE 和 p4vasp 批量打开所有的IMAGES中的CONTCAR,挨个查看。

小节

本节,我们主要是在前面计算的基础上,把gamma点的输出结果备份了一下,在此基础上,增加K点(提高精度)继续优化。最近瞎忙了很多事情,让很多人等急了,给大家拜个晚年,祝大伙猪年都学会麻溜滴算过渡态。

通过job-ID快速进入计算目录,可以有效地避免狂输cd命令,提高我们的计算效率。本节我们通过在国科智算超算中心的演示,介绍一个Slurm任务管理系统,保存任务ID,计算目录,并快速根据ID进入计算目录的方法。

核心思想:

有很多种方法可以实现:通过Job-ID进入计算目录的办法。这里只介绍本节的思路。

  • 通过squeue 命令获取用户所有的任务的ID;
  • 通过scontrol 命令获取所有ID对应的计算目录
  • 将计算ID,以及对应的目录保存到一个txt文件中。
  • 通过在~/.bashrc中使用scource命令,激活bash脚本中进入目录的操作。

之所以将计算目录保存到txt文件中,是因为如果计算已经完成的话,则不可以通过scontrol命令来获取目录,这时候,我们就可以在txt文件中查找相应的信息。

准备流程

按照流程操作,每一步都很重要,漏掉的话,可能会导致功能无法实现。

1 下载脚本:

ent.shjob_check.py复制到~/bin 目录下,赋予可执行权限。

1
2
gkzshpc101@login02:~/bin$ chmod  u+x ent.sh
gkzshpc101@login02:~/bin$ chmod u+x job_check.py

2 在bin目录下创建一个:job-check的目录,以及在该目录中创建一个空的job_list.txt文件,用来存储我们的任务信息。

1
gkzshpc101@login02:~/bin$ mkdir job_check && touch job_check/job_list.txt

3 修改job_check.py脚本中

1)13行中的账户信息,gkzshpc101改成你自己的账户

2) 32行中的那个目录信息,改成下面目录输出的。

1
2
3
gkzshpc101@login02:~/bin$ cd job_check 
gkzshpc101@login02:~/bin/job_check$ pwd
/public1/home/gkzshpc101/bin/job_check
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/env python
'''
This script is used to
1) record the job_ID and directories
2) print the job path from the job_ID
'''
from subprocess import Popen, PIPE
import numpy as np
import sys

def get_id():
list_j = [] # list for the job_ID
process = Popen(['squeue', '-lu', 'gkzshpc101'], stdout=PIPE, stderr=PIPE)
stdout, stderr = process.communicate()
list_out = stdout.split('\n')[2:]
for i in range(0, len(list_out)-1):
list_j.append(list_out[i].split()[0])
return list_j

def get_dir(job_id):
job_dir = None
command = 'scontrol show job ' + job_id
process1 = Popen(command, shell = True, stdout=PIPE, stderr=PIPE)
stdout1, stderr1 = process1.communicate()
for i in stdout1.split('\n'):
if 'WorkDir' in i:
#job_dir.append(i.split('=')[1])
job_dir = i.split('=')[1]
return job_dir

id_dict = {}
with open('/public1/home/gkzshpc101/bin/job_check/job_list.txt', 'a+') as file_in:
lines = file_in.readlines()
for line in lines:
key = line.split(':')[0]
value = line.split(':')[1]
id_dict[key] = value

list_j = get_id()

for i in list_j:
if i not in id_dict.keys():
id_dict[i] = get_dir(i)
file_in.write(i + ':' + get_dir(i) + '\n')

job_id = sys.argv[1]

for i in id_dict.keys():
if job_id in i:
print(id_dict.get(i).strip())

1
2
3
4
#!/usr/bin/env bash
pwd_work=$(job_check.py $1)
echo $pwd_work
cd $pwd_work

4 修改~/.bashrc文件,添加下面这一行:

1
alias ent='source ~/bin/ent.sh '

source 一下~/.bashrc文件: . ~/.bashrc

运行实例:

1
2
3
4
5
6
7
8
9
10
11
12
13
gkzshpc101@login02:~/test_ncore/test1/14$ ls
INCAR KPOINTS POSCAR POTCAR vasp.sh
gkzshpc101@login02:~/test_ncore/test1/14$ sbatch vasp.sh
Submitted batch job 14999
gkzshpc101@login02:~/test_ncore/test1/14$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
14999 operation vasp gkzshpc1 R 0:02 1 c0077
gkzshpc101@login02:~/test_ncore/test1/14$ cd
gkzshpc101@login02:~$ ent 999
/public1/home/gkzshpc101/test_ncore/test1/14
gkzshpc101@login02:~/test_ncore/test1/14$ ls
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT INCAR KPOINTS OSZICAR OUTCAR PCDAT POSCAR POTCAR REPORT vasp.out vasprun.xml vasp.sh WAVECAR XDATCAR
gkzshpc101@login02:~/test_ncore/test1/14$

国科智算提交VASP任务

今天QQ群里,有人问怎么在国科智算的超算中心提交VASP的任务。本着一言不合就写教程的态度,这一节我们就看下VASP的任务是怎么提交的。提交任务的脚本见群文件:vasp_qiangli.sh。下载后重命名成vasp.sh即可。想试用的,购买机时的可以加QQ群:608307988咨询一下。

Slurm

国科智算上系统采用的是slurm 任务调度工具。Simple Linux Utility for Resource Management,它是一个用于LinuxUnix 内核系统的免费、开源的任务调度工具,被世界范围内的超级计算机和计算机群广泛采用。它提供了三个关键功能。第一,为用户分配一定时间的专享或非专享的资源(计算机节点),以供用户执行工作。第二,它提供了一个框架,用于启动、执行、监测在节点上运行着的任务(通常是并行的任务,例如MPI,第三,为任务队列合理地分配资源。 大约60%的500强超级计算机上都运行着Slurm,包括2016年前世界上最快的计算机天河-2。 以上是来自维基百科的解释,具体的大家可以浏览Slurm 的官网:Slurm Workload Manager

Sbatch

在这个调度系统中,提交任务时我们需要用到命令: sbatch : https://slurm.schedmd.com/sbatch.html

sbatch 命令后面要跟一堆的参数,比如计算时间,节点数,邮箱,队列,调用的环境变量,任务名称等等。但这些信息通过命令直接输入又有些麻烦,所以我们把它们放到一个脚本里面,免得每次都重新输入一大长串的内容。而这个脚本,也就是我们本节内容的主角: vasp.sh

首先,我们浏览sbatch的详细参数: https://slurm.schedmd.com/sbatch.html

然后根据这些参数,我们就可以创建一个vasp.sh脚本了。

vasp.sh

下面就是vasp.sh 的主要内容:有2种写法,这两种写法也可以混着用,不影响。主要还是要参考sbatch的使用说明,需要什么就按照格式填写相应的内容。

写法(一)

1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/bash
#SBATCH -J BigBro ## Job Name
#SBATCH -o %j.out ## standard output
#SBATCH -e %j.err ## standard error
#SBATCH -p operation ## Partition
#SBATCH -N 1 ## Number of nodes
#SBATCH --ntasks-per-node=28 ## Each node has 28 tasks
#SBATCH -t 02-23:57:25 ## time for your job: 2 d,23 h ,57 min and 23 s

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std

写法(二)

1
2
3
4
5
6
7
8
9
10
11
#!/bin/bash
#SBATCH --job-name=BigBro ## Job Name
#SBATCH --output=%j.out ## standard output
#SBATCH --error=%j.err ## standard error
#SBATCH --partition=operation ## Partition
#SBATCH --nodes=1 ## Number of nodes
#SBATCH --ntasks-per-node=28 ## Each node has 28 tasks
#SBATCH --time=02-23:57:25 ## time for your job: 2 d,23 h ,57 min and 23 s

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std

上面每一行的含义,大师兄都注释出来了。

从第二行往下一次为:

  • 任务的名字,
  • VASP的标准输出,
  • 错误输出,
  • 任务运行的分区,
  • 使用的节点,
  • 每个节点的的tasks数目。这个tasks可以理解为每个节点的核数,国科智算的新机器每个节点是28个核。
  • 以及你的计算所用的时间。

最后两行为:我们调用的环境变量以及运行vasp程序。

需要注意的是,脚本里面的内容比如-N(修改任务所需的节点数目)、-J(修改任务的名字)这些我们频繁更换的,可以从脚本里面拿出来,在命令中运行。

1
sbtach -N 2 -J test vasp.sh 

本人自己的计算,一般28个核就够,任务名字也懒得修改。就把-N,-J写到vasp.sh中了。大家根据自己的任务特点自动修改就行了。下面我们具体演示一下。

实例操作1:

我们要在operation分区,运行一个VASP的单点计算,任务名称为:single, 使用2个节点,限制时间为2个小时;那么脚本的修改以及任务提交如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
gkzshpc101@login02:~/ex-A11/test_single$ 
gkzshpc101@login02:~/ex-A11/test_single$ ls
INCAR KPOINTS POSCAR POTCAR vasp.sh
gkzshpc101@login02:~/ex-A11/test_single$ cat vasp.sh
#!/bin/bash
#SBATCH -J single
#SBATCH -o out.%j
#SBATCH -e err.%j
#SBATCH -p operation
#SBATCH -N 2
#SBATCH --ntasks-per-node=28
#SBATCH -t 02:00:25

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std
gkzshpc101@login02:~/ex-A11/test_single$ sbatch vasp.sh
Submitted batch job 35525
gkzshpc101@login02:~/ex-A11/test_single$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
35525 operation single gkzshpc1 R 0:09 2 c[0032-0033]
gkzshpc101@login02:~/ex-A11/test_single$ ls
CHG CHGCAR CONTCAR DOSCAR EIGENVAL err.35525 IBZKPT INCAR KPOINTS OSZICAR out.35525 OUTCAR PCDAT POSCAR POTCAR REPORT vasprun.xml vasp.sh WAVECAR XDATCAR
gkzshpc101@login02:~/ex-A11/test_single$

实例操作2:

我们在Operation分区,运行一个VASPCI-NEB的过渡态计算任务,任务名称为:NEB,使用4个节点,限制时间为12个小时;那么脚本的修改以及任务提交如下:插了7个点,用112个核算,16个核算一个点。下面有点长,文字描述就到此结束,自己慢慢看,希望对大家有所帮助。注意:脚本里面,我们换成编译了VTSTvasp 5.4.1版本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
gkzshpc101@login02:~/ex-A11/test_neb$ ls
00 01 02 03 04 05 06 07 08 FS INCAR IS KPOINTS POTCAR vasp.sh
gkzshpc101@login02:~/ex-A11/test_neb$ cat vasp.sh
#!/bin/bash
#SBATCH -J NEB
#SBATCH -o out.%j
#SBATCH -e err.%j
#SBATCH -p operation
#SBATCH -N 4
#SBATCH --ntasks-per-node=28
#SBATCH -t 12:00:25

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/541_neb/vasp.5.4.1/bin/vasp_std
gkzshpc101@login02:~/ex-A11/test_neb$
gkzshpc101@login02:~/ex-A11/test_neb$ sbatch vasp.sh
Submitted batch job 35526
gkzshpc101@login02:~/ex-A11/test_neb$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
35526 operation NEB gkzshpc1 R 0:07 4 c[0056-0057,0060-0061]
gkzshpc101@login02:~/ex-A11/test_neb$ ls
00 01 02 03 04 05 06 07 08 err.35526 FS INCAR IS KPOINTS out.35526 POTCAR vasprun.xml vasp.sh
gkzshpc101@login02:~/ex-A11/test_neb$ ls *
err.35526 FS INCAR IS KPOINTS out.35526 POTCAR vasprun.xml vasp.sh

00:
POSCAR

01:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT WAVECAR XDATCAR

02:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

03:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

04:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

05:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

06:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

07:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

08:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR
gkzshpc101@login02:~/ex-A11/test_neb$
gkzshpc101@login02:~/ex-A11/test_neb$ head -n 10 out.35526
running on 112 total cores
each image running on 16 cores
distrk: each k-point on 16 cores, 1 groups
distr: one band on 1 cores, 16 groups
vasp.5.4.1 05Feb16 (build May 24 2018 19:36:47) complex

前面的准备工作完成了,INCARKPOINTSPOSCAR(IMAGES), POTCAR也检查完了。剩下的就是准备脚本提交文件了。本节主要是在天河II号超算中心上给大家简单示范一下:

  • 1)准备脚本,提交任务;
  • 2)过渡态任务运行时候的查看;
  • 3) VTST脚本nebresults.pl的安装和使用。

提交任务

首先通过命令操作,熟悉下天河II号超算中心提交NEB计算的流程。

1
2
3
4
5
6
7
8
9
10
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cat job_sub
#!/bin/bash
export LD_LIBRARY_PATH=/THFS/opt/intel/composer_xe_2013_sp1.3.174/mkl/lib/intel64:$LD_LIBRARY_PATH
yhrun -p gsc -n 24 /THFS/opt/vasp/5.4.4_neb/vasp.5.4.4/bin/vasp_std
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/$ yhbatch -p gsc -N 2 -J test job_sub
Submitted batch job 1004311
  • 每个节点有24核,我们使用了2个节点,共48核,来计算这个NEB任务。
  • 共有8个IMAGES,也就是6个核算一个IMAGE。
  • 我们调用的VASP版本是:/THFS/opt/vasp/5.4.4_neb/vasp.5.4.4/bin/vasp_std
  • 本例子下载链接:https://pan.baidu.com/s/1shjrJsRuDNYGEd2qDqvdFQ 提取码:k3pp

查看NEB计算

NEB一般收敛的很慢,需要跑很多离子步才会收敛,所以大家要心里有个准备,这也是为什么我们先用gamma点粗算一下,再用高密度的K点计算的原因,总之就是节约时间。虽然NEB要花很长时间才收敛,但提交任务后,我们不能守株待兔似的等着NEB算完,而是要勤检查结果,因为NEB的计算中,也会经常出现结构跑乱的情况。那么该怎么检查呢? 一看能量,二看结构,三看能量和结构。

一看能量:

主要是因为在terminal下面,相对于打开可视化软件查看结构来说,我们通过命令提取每个IMAGEOUTCAR能量信息的操作更加方便些而已。更重要的其实还是结构的变化,也就是你的反应路径。查看能量有2个办法:

自己用脑子想象或者写脚本:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ta.sh
01 -298.00705385
02 -297.96462657
03 -297.92754864
04 -297.90142637
05 -297.88475532
06 -297.87681445
07 -297.89025285
08 -297.93374021
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cat ~/bin/ta.sh
#!/usr/bin/env 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
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

用脑子想象:

ts.sh 脚本的输出中,从01的能量一直往下看,自己脑子里面自带一个xy的坐标系,将这些能量填到坐标系中,直至08。你会想象得到这样的一个曲线: 能量从01的-298.007慢慢上升到06的-297.87,然后再下降到08的-297.934。这个曲线的顶点再06处,也就是我们粗算得到的过渡态。

当然,你也写个小脚本提取每个IMAGE中的能量信息,然后画图即可。

使用VTST的脚本:nebresults.pl

前面Ex72中,我们讲过了怎么准备VTST的那些脚本,这一节,我们讲另外一个办法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
iciq-lq@ln3:/THFS/home/iciq-lq/bin$ cd vtst/
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ ls
vtstscripts.tgz
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ tar -zxvf vtstscripts.tgz
vtstscripts-937/
vtstscripts-937/pos2xyz.py
vtstscripts-937/sum_dos_np
vtstscripts-937/chg2cube.pl
vtstscripts-937/chgsplit.sh
vtstscripts-937/akmcprocess.pl
*
*
*
vtstscripts-937/nebbarrier.pl
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ ls
vtstscripts-937 vtstscripts.tgz
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ cd vtstscripts-937/
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ ls
2con.py center.py diffcon.pl dymanalyze.pl dymseldsp.pl insplot.pl ....
  • 使用vim打开nebresults.pl文件,将58-71行注释掉,或者可以使用sed命令:
1
sed -i '58,71s/^/#/g' nebresults.pl

效果如下图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
57 # Zip the OUTCARs again
58 #print 'Zipping the OUTCARs again ... ' ;
59 #$zip = $ENV{'VTST_ZIP'};
60 #if($zip eq ''){ $zip = 'gzip'; }
61 #
62 #$i = 0;
63 #$string = "00";
64 #while(chdir $string) {
65 # system "$zip OUTCAR";
66 # $i++;
67 # if($i < 10) { $string = "0$i"; }
68 # elsif($i < 100) { $string = "$i"; }
69 # chdir $dir;
70 #}
71 #print "done\n";
  • ~/.bashrc文件中设置脚本的目录:
1
2
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ pwd
/THFS/home/iciq-lq/bin/vtst/vtstscripts-937
  • 复制上面pwd命令的目录,打开~/.bashrc 文件,添加这一行:
1
export PATH=$PATH:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937
  • source一下~/.bashrc文件:
1
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ . ~/.bashrc
  • 进入计算的目录下,运行nebresults.pl命令,操作如下图:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

No OUTCAR in 00
Unziping the OUTCARs ...
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cp 01/OUTCAR 00
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

No OUTCAR in 09
Unziping the OUTCARs ...
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cp 08/OUTCAR 09
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.003505 -298.009300 0.000000
1 0.003505 -298.009300 0.000000
2 0.008701 -297.969500 0.039800
3 0.015508 -297.932900 0.076400
4 0.018300 -297.905700 0.103600
5 0.016765 -297.887500 0.121800
6 0.011627 -297.877000 0.132300
7 0.007293 -297.885300 0.124000
8 0.003607 -297.929100 0.080200
9 0.003607 -297.929100 0.080200

Extremum 1 found at image 0.210729 with energy: 0.002874
Extremum 2 found at image 0.788081 with energy: -0.002910
Extremum 3 found at image 6.265931 with energy: 0.132792
Extremum 4 found at image 8.211450 with energy: 0.074781
Extremum 5 found at image 8.788801 with energy: 0.085561

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

上面的那个图在Ubuntu的终端里面,通过gs mep.eps命令可以直接打开。

解释:

1) 使用nebresults.pl脚本的时候,我们需要在0009的文件夹中分别放上初始结构和末态结构所对应的OUTCAR。上面操作中,大师兄把0108中的OUTCAR分别复制到了0009中。

为什么这么做呢?

答:偷懒。原因有2个:

A)之前优化结构都是用到的高密度的K点,直接把OUTCAR拿过来用的话,会出现能量差别很大。做出图来能量很奇怪的样子。如下面的例子所示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
  iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/high_k$ ta.sh
00 -321.15187450
01 -298.00935617
02 -297.96955082
03 -297.93296493
04 -297.90579729
05 -297.88754704
06 -297.87700139
07 -297.88531984
08 -297.92917858
09 -321.09850728
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/high_k$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.023782 -321.151800 0.000000
1 0.003505 -298.009300 23.142500
2 0.008701 -297.969500 23.182300
3 0.015508 -297.932900 23.218900
4 0.018300 -297.905700 23.246100
5 0.016765 -297.887500 23.264300
6 0.011627 -297.877000 23.274800
7 0.007293 -297.885300 23.266500
8 0.003607 -297.929100 23.222700
9 0.011589 -321.098500 0.053300

Extremum 1 found at image 0.000009 with energy: -0.000000
Extremum 2 found at image 6.265997 with energy: 23.275310

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

看到了没有,上面例子中0108 结构的能量很大。这是因为参考的能量是00的。对于同一个体系,K点大的话,绝对能量更负一些。能量差别太大,导致过渡态的能量变化都可以忽略掉了。

B)当然我们也可以分别对0009对应对的机构算个单点。生成各自对应的OUTCAR,然后再使用nebresults.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ta.sh
00 -298.03233582
01 -298.00705385
02 -297.96462657
03 -297.92754864
04 -297.90142637
05 -297.88475532
06 -297.87681445
07 -297.89025285
08 -297.93374021
09 -297.97265721
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.100196 -298.032300 0.000000
1 0.003505 -298.009300 0.023000
2 0.008701 -297.969500 0.062800
3 0.015508 -297.932900 0.099400
4 0.018300 -297.905700 0.126600
5 0.016765 -297.887500 0.144800
6 0.011627 -297.877000 0.155300
7 0.007293 -297.885300 0.147000
8 0.003607 -297.929100 0.103200
9 0.084593 -297.972600 0.059700

Extremum 1 found at image 0.057295 with energy: -0.000155
Extremum 2 found at image 6.265997 with energy: 0.155771

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

仔细观察,你会发现,00010809 的能量差别基本不大。所以上面我们的偷懒做法也是可取的。

敲黑板:

一般在粗算的时候,不用考虑那么多细节,怎么粗怎么弄。所以,我们前面讲的这个懒办法是在粗算情况下的操作。由于高K点的计算我们在优化初末态结构的时候已经有了,所以提高精度的时候就可以直接用了。

本节要求

1) 会提交NEB的计算任务

2) 会使用自己写的或者nebresults.pl脚本查看能量信息。

下一节,我们学习怎么分析这些能量。

前面我们介绍了怎么准备过渡态计算的INCAR,POSCAR,POTCAR和POTCAR,这一节理所当然,就是要提交任务了。过渡态计算的时候,和大多数的计算是一样的:检查输入文件,准备好脚本,敲命令提交任务。所以,我们就在这三个方面展开,具体讨论一下里面的细节部分,避免不必要的出错,提高计算的成功率和效率。

心态

算过渡态的时候,要把自己想象成一个整装待命的狙击手,要熟悉自己的枪(VASP),周围的环境(影响计算的因素),以及狙杀的对象(计算的体系)。

也就是VASP的输入文件,开枪前,先检查一遍自己的枪有没有问题。

怎么检查呢?

大师兄本人的方法就是:不断在心里按着顺序默念:INCAR,KPOINTS,POSCAR,POTCAR。然后根据自己的经验检查并修改。这里你也许会说,大师兄,那是因为你有经验了啊,才能这么做,没有经验的我,该怎么办呢?

不要怕,我们下面慢慢分析。

INCAR

对于大多数的计算,影响计算时间的主要因素是KPOINTS,体系的大小,INCAR中关于收敛的参数。这三者紧密相连。而INCAR中的计算在粗算的这一个环节里面,对时间的影响不是很大。这给我们的启示就是:不管是粗算,还是后面提高精度的计算,我们使用同样的设置即可。大家一般来说的低高精度的计算,主要集中在EDIFF和EDIFFG这两个参数上。体系不是很特殊,EDIFF=1E-5; EDIFFG=-0.02 就足够了。完全没有必要再粗算的时候将它们的数值调大点或者调小点。如果可以这样思考的话,那我们就可以将INCAR的参数提前设置好,用的时候直接拿来就可以了。

想到了这一点,下面我们的任务就轻松了,可以写个脚本,直接敲个命令生成CI-NEB计算的INCAR。这里大师兄用的是另外一个办法:

1) 先准备好一个INCAR,里面肯定把CI-NEB计算的参数等乱七八糟的都写好了。

2)将INCAR放在一个固定的目录下面,比如:~/bin/INCAR_neb

3)提交任务前,将INCAR复制过来即可。

4)划重点:alias 的骚操作。

1
2
cp ~/bin/INCAR_neb INCAR 
alias neb='cp ~/bin/INCAR_neb INCAR'

这样的话,我们直接敲个命令,就可以得到INCAR文件了。

INCAR的检查:

前面的做法,并不代表我们就可以直接使用INCAR文件了。能不能直接用,取决于我们的具体计算细节。若INCAR中的参数与当前任务的不一致,或者缺少某些参数,我们就需要更新一下,主要有以下几个方面:

1) INCAR中默认的插点数目和实际不一致,记得更新。

2)若POSCAR中的结构有磁性,并且和INCAR中的不一致,记得更新;

3)如使用DFT-D2,INCAR中C6和R6的设置和POSCAR中的元素不一致,记得更新;

4)如使用DFT+U,INCAR中U和J的设置与POSCAR中的元素不一致,记得更新。

5)若计算的体系键能比较大,ISPRING这个参数可以稍微修改下,比如从-5调小到-8。

6)若发现slab体系粗算的时候,不怎么收敛,偶极校正取消一下。

KPOINTS的检查:

我们在粗算阶段,KPOINTS直接用gamma点即可。为了方便可以使用脚本kponits.sh 1 1 1 直接生成。也可以添加到我们的alias这个命令中:

1
alias neb='cp ~/bin/INCAR_neb INCAR && kpoints.sh 1 1 1 '

POTCAR 的检查:

可以使用我们前面的自动生成与POSCAR对应POTCAR的脚本:

1
cp 01/POSCAR . && pospot.sh && rm POSCAR 

为了方便,也可以直接将上面这句话添加到alias里面:

1
alias neb='cp ~/bin/INCAR_neb INCAR && kpoints.sh 1 1 1 && cp 01/POSCAR . && pospot.sh && rm POSCAR '

POSCAR的检查:

POSCAR被归结到枪这一部分,其实也是我们狙杀的对象。在粗算这部分的检查,需要注意的地方有:

1) 确保slab被固定住了;否则你算完会发现表面的原子都疯掉了。前面我们讲过怎么通过脚本或者sed命令批量固定了,就不再过多叙述。

2)使用p4v一次性打开所有的POSCAR文件:p4v 0*/POSCAR ; 这里大师兄也是用了alias

1
2
alias pp='p4v 0*/POSCAR'
alias ap='ase-gui 0*/POSCAR'

上面的2个操作都是在Linux下面的,p4v和ase-gui大家根据习惯随便选一款即可。

  • 使用p4vasp,然后挨个点IMAGE对应的POSCAR,检查结构对不对;

  • 使用ase-gui 将所有的文件全部打开后,摁键盘的PageUp和PageDown这2个键,类似于自己看一帧一帧的电影,体验也很棒。

  • 检查插点的时候,原子是不是乱套了,没有一一对应;
  • 检查我们插的点,是不是具有很好的化学或者物理意义,也就是预先评估下自己猜的这个路径靠不靠谱。

3) 使用Windows的用户先不用着急,大师兄正在学习VASPkit+VDM的骚操作。学会了下一节会及时更新,并提交过渡态的计算。

前面一节,我们提到说,在Ubuntu或者其他Linux系统下面,可以使用p4vasp或者ASE将结构批量打开,查看我们初步设置的NEB路径是否合理。但是在Windows系统下,我们不方便使用命令进行查看。这里介绍一下在Windows下面通过VASPkit结合VMD查看NEB路径结构的方法。

软件的获取:

vaspkit的使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$ ls
00 01 02 03 04 05 06 07 08 09 INCAR KPOINTS NEB.pdb POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$ rm NEB.pdb
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$ ls
00 01 02 03 04 05 06 07 08 09 INCAR KPOINTS POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$ vaspkit

+---------------------------------------------------------------+
| VASPKIT Version: 0.71 (16 Nov. 2018) |
| A Pre- and Post-Processing Program for VASP Code |
| Official Website: http://vaspkit.sourceforge.net |

*
*
*
*
0) Quit
------------>>
5
==================== Catalysis-ElectroChem Kit ==================
501) Thermal Corrections for Adsorbate
502) Thermal Corrections for Gas
503) d-Band Center (experimental)
504) Convert NEB-Path to PDB Format for Animation

0) Quit
9) Back
------------>>
504
+-------------------------- Warm Tips --------------------------+
See An Example in vaspkit/examples/neb_animation.
+---------------------------------------------------------------+
-->> (1) Reading Structural Parameters from 00/POSCAR File...
+---------------------------------------------------------------+
*
*
*
+---------------------------------------------------------------+
| Selective Dynamics is Activated! |
+---------------------------------------------------------------+
-->> (*) Reading Structural Parameters from 09/POSCAR File...
+---------------------------------------------------------------+
| Selective Dynamics is Activated! |
+---------------------------------------------------------------+
+---------------------------------------------------------------+
| * DISCLAIMER * |
| CHECK Your Results for Consistency if Necessary |
| Bug Reports and Suggestions for Improvements Are Welcome |
| Citation of VASPKIT Is Not Mandatory BUT Would Be Appreciated |
| (^.^) GOOD LUCK (^.^) |
+---------------------------------------------------------------+
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$ ls
00 01 02 03 04 05 06 07 08 09 INCAR KPOINTS NEB.pdb POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$
iciq-lq@ln3:/THFS/home/iciq-lq/demo_neb$

  • 在终端输入:vaspkit
  • 输入:5 后回车
  • 输入:504 后回车
  • 然后你会得到一个NEB.pdb文件。
  • pdb文件包含了00到09这几个文件夹中POSCAR的结构信息,用以VMD进行 查看。

VMD 查看pdb文件

使用VMD查看pdb的方法,我们将VASPkit中的具体说明拿过来展示一下:

  • Windows系统中启动VMD程序
  • 将我们在服务器中生成的NEB.pdb文件下载到本地,然后拖到VMD的界面
  • VMD主窗口选择菜单 Display —> Orthographic 正交显示模式
  • VMD主窗口选择菜单Graphics —> Representations —> Drawing Methods 选择 CPK
  • 默认是不显示盒子边界的,在VMD主窗口选择菜单 Extensions ,选择 Tk Console , 在弹出的VMD TkConsole 窗口中输入 pbc box -color white ,然后回车,查看模型结构。
  • 点界面的右下角的箭头后,你可以看到我们初步猜测的NEB路径中原子快速动起来了。箭头左面有个speed,我们可以调节原子的速度。
  • VMD主窗口选择菜单 Mouse —> Label —> 2, 然后去模型界面上,点与NEB路径中最相关的2个原子,就可以查看NEB路径中,原子间距离随着IMAGE结构的变化了。

总结

Windows下的用户,在做过渡态计算的时候,可视化是一个痛点,通常来说,Images中的结构都只能一个一个打开查看,计算的时候不能很好地体会一个反应的发生路径。使用VASPkit则可以顺利地解决这个问题,这个功能更详细的说明,请参考VASPkit的使用手册。当然,这个方法也适用于Linux操作系统下VASPkit + VMD的操作。此外,大师兄还是建议大家有余力的时候多多接触类似于Ubuntu,Centos这样的Linux操作系统。

过渡态计算INCAR的设置

vasp的过渡态任务所需要的五个文件:INCAR, KPOINTS, POSCAR, POTCAR以及脚本。我们已经学习了三个,还差INCAR和脚本,这一节,我们主要介绍INCAR中参数的设置。很多读者都没有接触过渡态的计算,总是幻想着过渡态是如何地高大上,令人可望而不可即。其实很简单,记住前面我们说的那句话,过渡态的计算无非就是优化一个特殊的结构而已。既然是优化结构,那么我们就可以在前面结构优化的INCAR基础上稍加修改,让它变成过渡态的INCAR即可。

修改一共有以下几个重要的部分:

第一部分

LCLIMB = .TRUE. 告诉vasp你要开始使用CI-NEB方法算过渡态了。

注意1: LCLIMB = .TRUE 这个写法是错的,因为TRUE后面少了个点!LCLIMB = T 这个写法是可以的。

第二部分

IMAGES = 8 告诉vasp,你插了8个点。VASP不是智能的,根据目录下文件夹的数目自动帮你数一下插点的个数,我们要自己设置一下。

注意2:INCAR中IMAGES的数目和实际的不符,这是很多人常犯的错误。比如插了8个点,INCAR中却是IMAGES= 4,这会导致VASP只读取00到05的结构,从而05到08文件夹中没有输出文件。

第三部分

过渡态计算的优化器(Optimizer)

CI-NEB是基于Force 也就是力来获取能量最低的反应路径。优化的方法有2个选择,一个是VASP默认的。一个是VTST中自带的。

3.1)VASP默认的优化方法:

IBRION =1 (quasi-Newton) 和3 (quick-min) 是基于力的优化方法。一般来说,如果你感觉自己插的点基本上就是反应的路径了,可以使用IBRION= 1。如果你的初始和末态结构不是很理想,插的点也是马马虎虎,那么3则是一个很好的选择。IBRION选择完了,我们还要设置一个合理的POTIM。个人经验0.1-1.0之间都是可以接受的。

所以:使用VASP自带的优化器:

IBRION = 1 或者3

POTIM = 取个合理的数值。

3.2)VTST中也自带了一些基于力的优化方法。

http://theory.cm.utexas.edu/vasp/optimizers.html#optimizers

使用VTST自带的优化方法,我们先要关闭VASP自带的,加下面2个参数。

1
2
IBRION = 3
POTIM = 0

然后你就可以通过IOPT来选择自己喜欢的了。下面红色部分是官网的介绍:

1
2
3
4
5
6
7
8
Optimizer input parameters
The following parameters are readfrom the INCAR file.
(IOPT = 0) Use VASP optimizersspecified from IBRION (default)
(IOPT = 1) LBFGS = Limited-memoryBroyden-Fletcher-Goldfarb-Shanno
(IOPT = 2) CG = Conjugate Gradient
(IOPT = 3) QM = Quick-Min
(IOPT = 4) SD = Steepest Descent
(IOPT =7) FIRE = Fast Inertial Relaxation Engine

注意:如果你已经打算用VTST自带的优化器,那么IOPT =0 不要选。

因为0对应的是VASP自带的。使用VTST的这个前提是我们把VASP的优化器给关掉了。但你关掉后,又设置IOPT=0,这是一个自相矛盾的选择,会导致计算过程中,你的原子纹丝不动,因为POTIM= 0.

所以,使用VTST自带的优化器:

IBRION = 3

POTIM = 0

IOPT = 1,2,3,4,7 有5个选择,官网建议使用1或者2。

第四部分

SPRING = -5 (这是默认值)

这个参数是干嘛的? 查阅一下VTST的官网,链接如下:

http://theory.cm.utexas.edu/vtsttools/neb.html

The nudged elastic band (NEB) is a method for finding saddle points and minimum energy paths between known reactants and products. The method works by optimizing a number of intermediate images along the reaction path. Each image finds the lowest energy possible while maintaining equal spacing to neighboring images. This constrained optimization is done by adding spring forces along the band between images and by projecting out the component of the force due to the potential perpendicular to the band.

什么意思呢?打个比方,这8个IMAGES就是一条绳上的8只蚂蚱,这些蚂蚱只能在一个方向上跳,在优化的时候,蚂蚱跳的太远,或者太偏就会被拉回来。拉回来的这个力就是通过SPRING这个参数来设置的。而我们之前的优化优化计算中,一条绳上只栓一只蚂蚱,该蚂蚱则比较自由,前后左右可以随便跳,并且没有人往回拽,这也是过渡态计算和普通优化所不同的地方。具体的理论部分,大家自行查阅NEB相关的参考文献。

如果你算的一个基元反应,两个原子之间的键很强,那么我们就需要将SPRING这个参数设置的更负一些,比如SPRING= -10,-15 或者-20。(不一定是-5的倍数,也可以是-6,-11 等。)如果你不知道怎么设置,一般来说默认值-5就足够用了。

第五部分

收敛相关的参数(粗算中的芝麻部分):

当然,粗算的时候,还要设置电子步以及以及离子步的收敛标准。

EDIFF = 1E-4 完全够用,不放心可以设置成-5;

NELM = 60 或者使用默认值40;

EDIFFG: 由于我们这一步是粗算,不用那么精确,EDIFFG = -0.05 完全够用了;

NSW: 个人经验,一般跑个50-60步左右就可以大体看个样子出来了,不用非得等计算达到EDIFFG设置的收敛标准。所以,前面的EDIFFG只是个形式罢了,这里我们先让这个参数占个坑,避免下一步提高精度计算的时候忘记这个参数。

小结一下:如果我们使用VASP自带的优化器,过渡态计算相关的参数如下:

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

如果你粗算了一下,跑了60步没有收敛,自己查看下各个IMAGES中的CONTCAR结构,感觉差不多就可以继续提高精度继续算了。如果感觉结构不是很好,那么就把各个IMAGES中的CONTCAR复制成POSCAR,然后继续再跑。后面大师兄会带着你一步一步跑,这个先不用着急。

总结:

这一节我们主要学习了过渡态计算中INCAR的一些基本设置。插点算过渡态时必须要有的参数如下:

1)告诉程序给我们要开始算过渡态了,并且指定插点的数目(这两个参数必须要写!!!)

LCLIMB = .TRUE.

IMAGES = 8

2)选择优化器(optimizer)

有VASP自带的和VTST自带的两种,大家选一种即可。

3)使用多大的劲控制蚂蚱不乱跳:

SPRING = -5 (默认值)

4)其他的参数跟前面结构优化的基本一致。

在准备工作中:第1)个是最关键的,不设置就不给你算。后面的2-4)由于都有默认的参数,即使不设置也不会导致VASP罢工的情况。但它们是一些经验性很强的参数,需要结合自己的体系具体设置。不过,经验性再强,目的也只有一个,就是过渡态要算准,所以一定要将能量和结构结合起来分析。尤其是结构,Linux用户多多使用类似

I) p4v 0*/CONTCAR

II) ase-gui 0*/CONTCAR

的命令将所有的结构一次性打开,分析结构变化和能量变化之间的关系。如果是在Windows下进行计算,自行搜索办法,或学习下一节中VASPkit + VMD的可视化办法。

前面2节我们讲了如何通过vtst的脚本进行插点,得到过渡态计算中特定数目的images,以及一个非常实用的小脚本来固定slab模型中底层的原子。本节我们就准备一下过渡态计算中:POSCAR、KPOINTS和POTCAR文件的准备工作。

关于Images

前面我们插了8个点,得到10个文件夹,分别为00, 01, 02,….09

细心的你会发现,每个文件夹中都会有一个POSCAR文件。00 中的POSCAR对应的是nebmake.pl IS FS 8命令中的IS,两个一模一样,同理,09 中的POSCAR是nebmake.pl IS FS 8命令中的FS。01到08位我们设想的一个反应路径。(通过p4vasp 0*/POSCAR 这个小技巧,来查看整个路径是不是合理)所有IMAGES中的元素行都是一样的,所以我们任取一个用来生成POTCAR文件。

结合之前VASP的使用经验,每个文件夹中都有一个POSCAR,那么算NEB的时候:是不是各个文件夹也中都要有INCAR, KPOINTS和 POTCAR以及提交任务的脚本呢?

如果你这么想,是很正常的。但我们算过渡态的时候,却不是这样操作的。

实际上,我们将00 到 09 这10个文件夹作为一个整体,看作一个POSCAR。和在这些文件夹相同的目录下, 有INCAR, KPOINTS 和 POTCAR文件。

如果你还不理解,看下面的两个图:首先我们放一个错误的准备工作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/wrong$ ls
00 01 02 03 04 05 06 07 08 09
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/wrong$ ls *
00:
INCAR job_sub KPOINTS POSCAR POTCAR

01:
INCAR job_sub KPOINTS POSCAR POTCAR

02:
INCAR job_sub KPOINTS POSCAR POTCAR

03:
INCAR job_sub KPOINTS POSCAR POTCAR

04:
INCAR job_sub KPOINTS POSCAR POTCAR

05:
INCAR job_sub KPOINTS POSCAR POTCAR

06:
INCAR job_sub KPOINTS POSCAR POTCAR

07:
INCAR job_sub KPOINTS POSCAR POTCAR

08:
INCAR job_sub KPOINTS POSCAR POTCAR

09:
INCAR job_sub KPOINTS POSCAR POTCAR

这就是刚刚说的,把INCAR, KPOINTS, POTCAR 和脚本都复制到images对应的文件夹中。知道了错误的准备文件,那下面我们看一下正确的做法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/wrong$ cd ../right/
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/right$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/right$ ls *
INCAR job_sub KPOINTS POTCAR

00:
POSCAR

01:
POSCAR

02:
POSCAR

03:
POSCAR

04:
POSCAR

05:
POSCAR

06:
POSCAR

07:
POSCAR

08:
POSCAR

09:
POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex74/right$

把0{0..9}这些文件夹作为一个整体,看作POSCAR。这些文件夹 和INCAR,KPOINTS, POTCAR,以及提交任务的脚本放在同一个目录下。

备注:

1) 上面本人插了8个点,纯属自己的日常计算习惯,并不属于什么过渡态的技巧。大家根据自己服务器的运算能力,设置合理的数目。

2) POTCAR与这些images中的元素数目要对应,这个大家应该都会了。

3) KPOINTS: 划重点

如果你从网上搜索相关的资料,很多时候会发现大家经常说,先粗算一下(或者预优化一下),然后再逐渐提高精度。算过渡态的时候也是这样子的,我们不能100%确认自己预想的反应路径(Images文件对应的结构)就是合理的,或者刚刚插的这些点结构都很粗糙,直接用高精度的计算会造成资源的大量浪费,也就是费力不讨好。再严肃点,纳税人的钱不是让你这样随随便便烧的。所以,我们要掌握一些合理的技巧,争取用最少的计算资源,最短的时间,来获取最合理的过渡态结构。前面说到了粗算,到底什么是粗算?粗算需要改什么参数?

这里我们需要了解芝麻和西瓜的故事

通过前面的学习,我们知道KPOINTS以及slab(计算体系)的大小对于计算的影响最大,它们就是西瓜。

而其他的参数,比如收敛标准,阶段能的大小等,都是小芝麻(当然,你收敛精度设置地极其高那也是个别情况,对于粗算的时候我们不考虑)。

说白了, 就是要抓住影响机时的关键因素。所以粗算的时候,我们需要改的主要有2个地方: KPOINTS 和 slab的模型(集中力量啃西瓜)。

Slab的修改难度较大, 对于新手不太实用,老司机也经常翻车,故先不做讨论。这里我们主要介绍通过调节KPOINTS的一个技巧来实现快速粗算的目的。

先声明:

以下适用于slab表面上的反应计算。如果你算的是体相里面的,比如原子的扩散等,只可作为参考,不能照搬。

在粗算的时候,建议使用低密度的的K点。上来就用高密度的K点计算,即使是土豪,也是不建议的。那么K点密度怎么个低法呢? 直接来个极端就好了,用gamma点。前面我们还讲过,如果用gamma点计算的时候,需要把slab的原子固定住。

为什么呢?

你可以使用4层的p(2x2)的Cu(111)的slab模型,放开上面2层,然后gamma点优化一个CO在上面的吸附,最后会发现表层的Cu原子驰豫得如同波浪一般起伏。也就是用gamma点会导致表层原子过分地驰豫,从而造成错误的结算结果。

但这也不是绝对的,如果你的slab在xy方向的尺寸很大,使用gamma点的时候满足我们之前讲过的 k*a 的经验规则,即使不固定原子,也是很安全的。(为保险,需要自己测试)如果slab在xy方向的尺寸很小,使用gamma点的就一定需要把原子固定住。固定slab的时候,这里面就有很多窍门了。

A)首先给大家介绍一个最简单的办法。也就是上一节的内容,使用塔马斯的POSCARtoolkit.py脚本,固定和放开任意slab里面的底层原子。比如,目前我们的例子,4层的Ru原子上面,H原子的扩散。我们就可以通过一个for循环来实现slab原子的快速固定。为了方便大家,写了一个基于POSCARtoolkit.py 的bash脚本,该脚本命名为fix.sh, 内容如下:

1
2
3
4
5
6
7
8
9
10
11
12
#!/usr/bin/env bash
## To use it
## fix.sh 4 POSCAR will fix the bottom 4 layers of POSCAR in all image filders
## fix.sh 2 CONTCAR will fix the bottom 2 layers of CONTCAR in all image folders
## The old POSCAR or CONTCAR will be replaced.

for i in 0* ; do
cd $i;
echo "$1" | POSCARtoolkit.py -i "$2" -f;
mv ${2}_C "$2"
cd $OLDPWD
done

注意:

1) 这里大师兄将默认的阈值设置为0.5 Å(修改的脚本)。因为原来1.5 Å有点大, 使用的时候会将表面的H和表层的Ru看作一层。

运行到这里的时候, 我们将所有Images中POSCAR的slab对应的原子全部固定住了。也就是粗算的POSCAR设置好了。

2) 不要偷懒直接复制代码就运行,由于网页上的格式问题,可能直接复制过来运行的时候回出错。

B) 除了使用脚本外,还有另外一个方便的办法

在课题刚刚开始的时候,就把slab模型中的原子按照层数(也就是z方向的大小)进行排列,前面我们讲过怎么sort坐标了,这里就不再啰嗦。

排列后的一个优势就是各个层的原子在POSCAR中的序号是连续的,可以很方便地进行选择然后通过使用sed 命令将POSCAR里面的T T T 批量转化为 F F F.

比如本例中:10-45为Ru原子,10-18, 19-27, 28-36, 37-45 从上往下数,分别对应第一、二、三、四层的Ru原子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
Ru   H
1.00000000000000
8.1377999784000004 0.0000000000000000 0.0000000000000000
4.0688999892000002 7.0475415119999996 0.0000000000000000
0.0000000000000000 0.0000000000000000 21.5631999968999999
Ru H
36 1
Selective dynamics
Direct
0.2212968850718230 0.2239409561475650 0.2936169649913790 F F F
0.2210036114118170 0.5561663540453879 0.2936884876706930 F F F
0.2212954084343650 0.8880992570872690 0.2936171801115680 F F F
0.5547059120221930 0.2208073369214740 0.2945957656788490 F F F
0.5547051050620690 0.5578229108467880 0.2945953069057910 F F F
0.5572530259134471 0.8880427972121010 0.2935860240494040 F F F
0.8910551675972660 0.2211397177458350 0.2944546614649750 F F F
0.8893412673473050 0.5559102687354450 0.2935156586054220 F F F
0.8893404937786400 0.8880862995751930 0.2935168781920150 F F F
0.1119432861914480 0.1106972088633780 0.1977748273385700 F F F
0.1112675774095350 0.4441060896203480 0.1977506775579700 F F F
0.1112675246777360 0.7779649556723670 0.1977514668340350 F F F
0.4440066349954530 0.1106943439331120 0.1977106006833940 F F F
0.4446549772136020 0.4443417041956370 0.1981365664070010 F F F
0.4440061461204080 0.7786377563630820 0.1977110173547200 F F F
0.7777009262273890 0.1113226147079440 0.1980782640541890 F F F
0.7777006933574599 0.4443146361556100 0.1980780083860340 F F F
0.7774653980933121 0.7779366397494840 0.1977260420610920 F F F
0.2222233413377670 0.2222233413358840 0.0992873228095945 F F F
0.2222233413377670 0.5555566746692190 0.0992873228095945 F F F
0.2222233413377670 0.8888900080025550 0.0992873228095945 F F F
0.5555566746710950 0.2222233413358840 0.0992873228095945 F F F
0.5555566746710950 0.5555566746692190 0.0992873228095945 F F F
0.5555566746710950 0.8888900080025550 0.0992873228095945 F F F
0.8888900080044310 0.2222233413358840 0.0992873228095945 F F F
0.8888900080044310 0.5555566746692190 0.0992873228095945 F F F
0.8888900080044310 0.8888900080025550 0.0992873228095945 F F F
0.1111100016643290 0.1111100016603930 0.0000000000000000 F F F
0.1111100016643290 0.4444433349937280 0.0000000000000000 F F F
0.1111100016643290 0.7777766683270571 0.0000000000000000 F F F
0.4444433349976650 0.1111100016603930 0.0000000000000000 F F F
0.4444433349976650 0.4444433349937280 0.0000000000000000 F F F
0.4444433349976650 0.7777766683270571 0.0000000000000000 F F F
0.7777766683310010 0.1111100016603930 0.0000000000000000 F F F
0.7777766683310010 0.4444433349937280 0.0000000000000000 F F F
0.7777766683310010 0.7777766683270571 0.0000000000000000 F F F
0.6419725686006360 0.3456803556224560 0.3432239557346220 T T T

在上面POSCAR基础上,放开上面的2层:

1
sed –i  '10,27s/F/T/g’ 0*/POSCAR

固定所有的Ru原子,可以:

1
sed –i  ’10,45s/T/F/g’ 0*/POSCAR

大家随意练习折腾,本人也相信,如果从本书的Ex-0一直认真学下来,如果sed 命令还不会的话,说明你真的没有认真在学习。本节的例子以及脚本下载链接:

https://pan.baidu.com/s/1qqySBpqc0zMu1lvpLVZuWA 提取码: 3ryx

总结

本节,我们介绍了准备文件中POSCAR, KPOINTS,和POTCAR,以及他们在目录里面的排列方式。为减轻计算压力,过渡态可以先粗算一下, 通过使用gamma点结合固定slab的方式来减少计算量。通过POSCARtooklit.py或者sed命令实现批量固定slab的方法。

在学习下面neb的计算前,思考了很久,感觉有必要先给大家介绍一个非常实用的python小脚本,经历过反复的修改,最终在群友小塔玛斯的努力下,完成了一个稳定的过渡态计算中非常实用小脚本:POSCARtoolkit.py

本节我们介绍一下这个脚本的具体使用方法。

版本和版权问题:

版本:

这个脚本有2个版本,分别适用于python3和python2的2.6及以上版本。

大家根据自己使用的python版本进行下载,如果运行的时候出现错误,请先检查是不是版本的问题。链接:https://pan.baidu.com/s/1JJHVM27NdIC-uZlxzlg-Mg 提取码:mt7e

版权问题:

这个脚本大家随便实用,如果感觉不错,碰到作者时,请他吃顿饭。作者坐高铁买不起盒饭时,在任何QQ群里呼救时,使用者需慷慨解囊,集资把盒饭钱凑够。

如果不接受这个条款,请自觉不使用该脚本。

使用前的准备工作

2.1)下载脚本,将脚本命名为:POSCARtoolkit.py

2.2)可执行化:

1
chmod u+x POSCARtoolkit.py

2.3)将脚本移到~/bin 文件夹下面:

1
mv POSCARtoolkit.py ~/bin

注释:

A) 2.2 和2.3步的顺序可以颠倒,不影响使用;

B) 如果没有~/bin文件夹,那么 mkdir ~/bin 手动创建一个即可。

坐标转换

这个脚本可以实现分数坐标(Direct)向笛卡尔坐标(Cartesian)的转换。VASP的输出结果是以分数坐标的形式存在CONTCAR中,但我们在操作模型的过程中,移动原子都是以Å来进行的,所以将分数坐标转化为笛卡尔坐标对于搭建模型很有帮助。脚本用法:POSCARtoolkit.py -i POSCAR

注意:

如果你的POSCARtoolkit.py脚本在和POSCAR一个目录下,使用下图中的命令,如果你已经完成了前面2中的步骤,脚本前面的python 不用输入。

描述1) -i 代表 input的意思, 后面紧跟你要转化的文件,可以是POSCAR,CONTCAR,也可以是其他的VASP的坐标文件。比如你把POSCAR命名成bigbro, 也是可以直接转换的。 POSCARtoolkit.py –i bigbro 即可。

描述2):脚本将含有分数坐标的POSCAR 转换成笛卡尔坐标的POSCAR_C。输出的文件被命名为POSCAR_C。

描述3)如果你预先把POSCAR命名为bigbro, 那么运行程序后,输出结果为 bigbro_C.

描述4)如果POSCAR已经为笛卡尔坐标,则转换终止,大家可以尝试一下转换一个笛卡尔坐标的POSCAR。

原子层数的固定

除了可以实现坐标系的转换,该脚本还可以根据Z方向的原子坐标,固定底部几层原子。

这在slab模型的相关计算中非常实用,如下图所示,在slab模型中,我们经常将底部的原子层固定,只放开表面的来进行计算。

比如我们把底部的三层原子进行固定。脚本的用法如下:POSCARtoolkit.py -i POSCAR_C -f

描述1) 脚本会根据阈值(默认1.5 Å)划分层,这里的1.5 Å指的是层间距。

描述2) 用户可以通过在脚本后面增加参数 -y 1.0 自己定义更小的阈值

描述3) 如果不想每次使用的时候使用参数-y,用户也可以直接在脚本里面修改阈值的大小,如下如,在第40行中,默认的阈值1.5被修改0.5。

当然,修改脚本默认的参数后,如果使用中使用的 –y, 还是会按照 –y 后面的阈值进行操作。

描述4) 输入完命令,回车后,会提醒用户输入固定的层数,这里我们输入3,回车,即固定底下三层原子,而其他原子放开。

描述5) 上图命令中,我们对POSCAR_C 进行操作,输出文件问POSCAR_C_C。

固定和放开用户选择的原子

除了固定层数外,该脚本允许用户选择部分原子放开,固定或者部分放开。这个功能实现的前提是POSCAR 或者 CONTCAR 中有 Selective dynamics 信息。如果没有,则可先固定任意的原子层(POSCARtoolkit.py -i CONTCAR -f ),这样的话Selective dynamics就会被写入到输出文件POSCAR_C中。然后用户再对POSCAR_C进行原子选择性操作。 下面我们先详细介绍一下脚本的用法,然后再加一些实例的操作来帮助用户理解。

用法:POSCARtoolkit.py -i POSCAR_C [-f or -r] -s [your selections]

描述1):POSCARtoolkit.py -i POSCAR_C 我们对输入文件POSCAR_C 进行操作

描述2): -f 和 –r 配合后面的 –s 进行操作。

  • -f 表示表示固定(fix)选中的原子,
  • -r 代表放开(relax)选中的原子,由于原子在xyz三个方向上都可以选择放开,所以使用-r的时候要配合F T 来进行操作。 如下:
    • -r FFT 代表只放开z方向,同理 –r TTF, -r TFT, -r TTT, -r FTT这些你就知道是怎么回事了。
    • -r 后面的FFT这三个字母之间可以有空格,也可以没有。也就是说:FFT 和 F FT, FF T 以及 F F T 效果是一样的。

描述3) -s 选项表示选择部分原子, 后面是你要选择的原子,选择项如下:

  • all 表示选中所有原子

  • 1-5 6 9 表示选中 第1-5个和6,9 号原子

  • Pt 表示选中所有的Pt原子

  • 1-5 6 9 Pt 表示选中 第1-5个和6,9 号原子和所有的Pt原子

描述4)其他未选择的原子,限制信息保持不变。

实例操作1: 如果我们想固定POSCAR_C中所有的原子:

1
POSCARtoolkit.py -i POSCAR_C  -f  -s  all

注释: -f 表示fix, -s all 表示选择所有的原子

POSCAR_C 中必须有Selective dynamics这一行

实例操作2: 将所有的原子只在z方向上放开:

1
POSCARtoolkit.py -i POSCAR_C  -r FF T  -s  all

注释:放开z方向上使用 –r FF T, FFT之间有无空格均可。

实例操作3: 将所有的C,H, O 原子在xy方向上放开:

1
POSCARtoolkit.py -i POSCAR_C  -r TTF  -s  C H O

实例操作4: 将所有的Pt原子和40号原子在z方向放开:

1
POSCARtoolkit.py -i POSCAR_C  -r  FFT  -s  Pt 40

批量转化

前面介绍的都是针对一个文件进行操作,由于固定原子层数的功能我们设置了一个交互,需要用户指定需要固定的层数(z轴从下到上)。

但这对于层数批量固定操作来说是个累赘,我们不想每操作一个文件就输入一次层数。

所以,如果想跳过交互,可以通过管道连接符 | 实现。比如固定底部4层:

1
echo 4 | python POSCARtoolkit.py -i CONTCAR –f

知道了这一点,我们就可以通过一个for循环进行批量固定层数的操作了。比如我们有POSCAR1, POSCAR2 到POSCAR100个文件,

1) 我们想批量固定它们的底部3层。

1
for i in POSCAR{1..10}; do echo 3 | python POSCARtoolkit.py  –i  $i  –f ; done 

2) 我们想批量将C H O原子在z方向上放开,zy方向上固定。(这个不需要管道连接符)

1
for i in POSCAR{1..10}; do python POSCARtoolkit.py  –i  $i  –r FFT –s C H O ; done

评价

这个脚本在模型操作的过程中,实用性很强,也很方便。实在是VASP计算中的一大利器。本脚本是群友tamas-zju-VASP(小塔玛斯)费煞苦心完成的。如果感觉不错,欢迎打赏。下载链接: https://pan.baidu.com/s/1JJHVM27NdIC-uZlxzlg-Mg 提取码:mt7e

不好意思,让大家久等了。最近一直很忙,后面还是会很忙。为了不让大家失望,抽时间把简单的过渡态计算讲解一下,这一节我们需要学习的是通过CI-NEB(Climbing Image Nudged Elastic Band)计算H原子在Ru(0001)表面上的扩散过程。

NEB计算过渡态的准备工作

准备工作1:

说到NEB或者CINEB,不得不提的就是Henkelman课题组的VTST:http://henkelmanlab.org/ . 我们需要下载VTST的code,然后将VASP编译一下,得到一个可以使用vtst的VASP版本。

链接:http://theory.cm.utexas.edu/vtsttools/download.html

很不幸,怎么编译,我一窍不通。如果这一关过不了的话,后面的也只能看看热闹,不能亲自实践了。如果你已经编译好了VTST的VASP版本, 那么就可以继续下面的学习了。

准备工作2:

我们将VTST的一些实用小脚本下载下来:http://theory.cm.utexas.edu/vtsttools/download.html

下载VTST Scripts后解压,把所有脚本都复制到 ~/bin 目录下:效果如下:

很多脚本估计你一辈子都不会用到,没关系,它们很小,不占存储空间,静悄悄躺在bin文件里面就行了。Linux老手们嫌烦可以把这些脚本放到一个文件夹里面,然后bashrc文件里面修改下路径就可以了。

准备工作3

修改nebresults.pl 文件,将57到71行注释掉,注释掉就是在每一行的开头加个#号,效果如下:

为什么需要这么做?

如果neb运行的时候,你使用这个命令,它就会把你的OUTCAR文件压缩,VASP找不到OUTCAR,就不知道该往哪里存储,然后就挂掉了。大家可以在neb任务运行的时候分别运行下注释前后的脚本,对比下就清楚了。

扩散的基础知识

准备工作完成,我们加深下扩散(diffusion)的一些概念:

2.1 扩散我们高中的时候应该就学习过了。具体的定义这里不讲了。大家自行参考教科书或者维基百科:链接如下:

https://en.wikipedia.org/wiki/Diffusion

对于分子在表面上的扩散,就跟我们在校园里瞎晃悠是一样的。哪天指不定遇到个对上眼的美女,过渡态就要开始了。

2.2 具体到分子在金属表面上的扩算,一般来说,扩散的能垒差不多是吸附能的12%左右。这里我们说的吸附能,指的是最稳定的那个结构所对应的。

为什么是12% 呢? 这只是一个经验性的结论,大家可以参考大牛Manos Mavrikakis的俺狗娃文章: (不要留言问我要这个文章…)

A Simple Rule of Thumb for Diffusion on Transition-Metal Surfaces

https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.200602223

2.3 从结构角度来看,举个面心立方的金属(111)面为例,表面上有fcc, hcp, top, 和 bridge这四个位点。如果fcc或者hcp的位点对应稳定的吸附构型,那么bridge的这个位点就是扩散的过渡态。所以,如果你想简单算一个扩散的过渡态,直接算bridge的吸附就可以了。但是,bridge的位点也不是那么容易就可以算出来的。前面我们讲过Cu(111)表面的吸附计算,很多人就反应bri的结构优化不出来。Bridge 的位点就跟下面的这个独木桥类似,少一不小心就掉下去了。

换个角度来分析:既然我们直接优化得不到bridge的结构,而这个结构恰恰就是过渡态,那么我们就可以使用算过渡态的方法来得到bridge的结构。当然,能直接优化出来bridge的结构是最好的。

过渡态计算的步骤

3.1 第一步:优化初始和末态结构。

首先我们先优化一个扩散前后的两个结构: H在FCC和HCP位点上的吸附。这个对于认真练习过前面计算的筒子们来说就是小菜一碟了。直接列出来top view 的示意图:

FCC site

HCP site

3.2 第二步:准备NEB的结构

这里我们需要用到VTST官网的一个小脚本:nebmake.pl,使用方法如下:

nebmake.pl IS FS N

1) 敲命令nebmake.pl

2) IS 指的是初始结构:

3) FS指的是末态结构

4) N 指的是你要插点的个数。

细节1:

IS, FS 是VASP的POSCAR或者CONTCAR。

可以是其他目录里面的POSCAR或者CONTCAR, 也可以是当前目录下的POSCAR或者CONTCAR。

你可以把初始结构的POSCAR命名成你的名字:bigbro , 末态命名成:bigbra. 运行的时候命令应该这么敲:nebmake.pl bigbro bigbra 8

下面图里面的三个事例,是等价的。自己随意领会:

细节2:

插入的点数要保证可以被使用的核数整除。比如我们打算用24个核进行计算,那么N可以是下面的几种情况:

N

Cores used for each Image

Total

1

24

24

4

6

24

6

4

24

2

12

24

12

2

24

3

8

24

8

3

24

注意:

表中最后一个很可能会出错,因为你用3个核来计算一个image. 上面说的不是绝对的,具体要根据自己的服务器和习惯来设置。

本人使用4个节点(每个节点12个核,共48个核)进行计算,一般N设置位4或者8.

细节三:

同样核数下,N设置的越大,计算每个image所需要的核数就会减少,导致计算变慢。

N需要怎么设置才好呢?

一方面需要大家去咨询一些现成的经验:问问师兄师姐,没有师兄师姐,就在现有的条件下自己去摸索摸索。当然这需要耗费大量的时间和精力,也不符合当下很多人急于速成的心态。

我们着重讲解另一个方面,也就是你对自己研究体系的掌控性。个人感觉分两点:

1)反应过程结构变化的理解,这需要我们在前面初末态的优化上下功夫;

2)化学键的理解:对一个化学键来说,它断裂时候的过渡态键长应该在这个化学键的1-2倍之间。1 指的是这个化学键本身,2指的是这个键被拉长了2倍。你肯定会说我在扯蛋,如果更经验性一点,应该在1.5倍左右。掌握这个这有助于你第一眼去判断自己的过渡态对不对。当然,1.5只是个粗糙的数值,不同的键会有不同的经验性数值。

3.3 第三步: 检查初末态结构的原子坐标是否是一一对应。

这是很多人算过渡态经常忽略的一步,有时候也是很费精力的一步。一个人是不是在闭着眼瞎算,从这一步基本上就可以看出来了。这里大师兄推荐一个linux下检查结构的方法。 命令:p4v 0*/POSCAR 一次性打开所有的Image结构,然后逐个点开,查看整个反应轨迹进行检查。

扩展练习

本节我们需要做的很简单

1) 浏览VTST的网站,阅读相关过渡态计算的步骤以及CI-NEB相关的文献;

2) 下载VTST的脚本,复制到~/bin文件夹,修改nebresults.pl脚本;

3) 计算H在FCC 和HCP的吸附;

4) 使用脚本生成NEB计算的Images文件。

5) 使用脚本生成前面两节:NH3翻转,以及乙烷旋转的Images文件。

总结

本节应该够新手们练习一阵子的了,下一节,我们介绍怎么把NEB的计算运行起来。

前面一节我们讲了直接优化NH3翻转的过渡态。这一节我们继续直接优化过渡态的介绍。计算乙烷分子(CH3CH3)绕C—C键旋转的过渡态。我们知道:乙烷有交错式和重叠式构型。而重叠式能量较高。想象乙烷分子绕着C—C旋转,从交错式到重叠式,再到交错式这一个过程。很容易就得到这个结论:重叠式的结构就是过渡态。所以,如果计算这个过渡态的工作就转化成了优化重叠式的构型。这也是本节的一个主要思想。


手动搭建初始的乙烷构型。

首先,我们要搭建一个初始的乙烷构型。把C—C键的两个C平行于z轴。思考一下:

1) 我们为什么要这样做?

2) 乙烷构型网上到处都是,随便下载一个就可以直接拿来使用,如果将任意取向的结构转化为C—C 键平行于z轴的结构呢?

下面是我们搭建的乙烷结构,并优化的结果。优化过程中,两个C原子在xy方向的坐标被固定住了。

很显然,我们的初始结构经过优化后转变成了交错式的构型。


搭建过渡态的初始构型

在重叠式的结构中,乙烷分子的两个CH3关于穿过CC轴的xy平面是对称的。所以,我们将上面CH3中H原子的xy坐标修改的和下面部分一样。

注意:

1) 只修改xy坐标,z坐标保持不变

2) 观察图中:3和8, 4和7, 以及5 和6 的xy坐标。

3) 固定xy坐标,直接优化z方向

下面是优化完的结果:计算很快就收敛了。


对比一下交错式和重叠式C—C的键长:

在重叠式中,C—C 稍微拉长了一点。

频率计算

在频率计算的时候,本人把所有的原子都放开了,如下图:

通过命令发现,有个256 cm-1的虚频。使用Jmol可视化一下:

过渡态计算完毕。


总结

大家等了很久,看完肯定会对本节失望。但过渡态计算里面的一些技巧很多都是自己琢磨出来的。比如我们将分子沿着平行z轴的方向放置,根据对称性手动搭建过渡态的结构,以及固定坐标,选择性地优化等。还是那句话:过渡态只是个结构,我们需要做的就是通过各种各样的办法来实现快速优化这个结构的目的。

本节计算例子下载链接: https://pan.baidu.com/s/18TdaLhWnz4_IDMjJJGBMFw 密码: ugjt

从本节起,我们开始学习下过渡态相关的计算。在进行下面的讲解之前,大家务必要记住这两点:

  • 1) 过渡态只是一个结构而已,过渡态的计算并不是你想象的那么神秘;也就是计算很容易。

  • 2)物理化学,结构化学书中的过渡态相关的基本知识要掌握,也就是你的理论基础要扎实。感觉基础不好的,认真看书。

过渡态计算的例子

我们遵循一个从简单到复杂的顺序,慢慢介绍一下表面化学中过渡态的一般计算思路和方法。简单的例子,本节我们选取的是VASP官网的经典计算:NH3的翻转,如下图:

如果你不太懂这个过程,可以把NH3看做一把雨伞。翻转的过程,就相当于大风把伞吹翻。

对于这个计算的一些参考资料,大家可以点下面这两个链接:

1) http://www.vasp.at/vasp-workshop/tutorials/tutorial_ammonia_flipping.pdf

Vasp官网的一个PPT介绍,这个PPT有些年头了,不过对于我们学习还是很有帮助的。建议看完PPT里面的内容(最好是自己练习一遍后),再继续下面的学习。

2) https://cms.mpi.univie.ac.at/wiki/index.php/Category:Examples

VASP wiki网页中关于过渡态计算的例子。


过渡态计算:初末态结构的优化。

过渡态计算的第一步:我们要知道自己想算什么,需要准备什么?

算过渡态,肯定指的是某一个反应或者过程的过渡态。过渡态的两边也分别对应着反应物和反应产物。所以,算一个反应的过渡态,不管计算过程如何,我们最终都会有三个结构:初始和末态结构,以及过渡态。初末态的结构和前我们已经讲过的优化过程是一样的,我们先回顾一番。


初态结构的优化

1)搭建结构(POSCAR)

该图中的三个结构分别对应的是NH3翻转的初始,过渡态以及末态结构。这里反应物和产物都是NH3。首先,我们肯定要有一个NH3的结构。这里大师兄就直接复制PPT里面的坐标快速制备一个POSCAR。

注意: PPT里面没有写元素符号,大家记得自己加上。下图是本人复制后制作的POSCAR。 (添加了元素行,Cartesian坐标。)

1
2
3
4
5
6
7
8
9
10
11
12
13
ammonia flipping
1.0
+9.0000000000 +0.0000000000 +0.0000000000
+0.0000000000 +10.0000000000 +0.0000000000
+0.0000000000 +0.0000000000 +11.0000000000
H N
3 1
Selective
Cartesian
+3.8185740000 +3.9721220000 +4.3936400000 T T T
+3.0000000000 +2.5542720000 +4.3936400000 T T T
+2.1814260000 +3.9721220000 +4.3936400000 T T T
+3.0000000000 +3.5000000000 +4.0000000000 F F F


2) 准备输入文件(INCAR)

优化NH$_3$的初始结构:(回顾气相分子计算需要注意的地方)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
System = NH3
ISMEAR = 0
SIMGA = 0.01

ALGO = FAST
ENCUT = 450
EDIFF = 1E-5

NWRITE = 0
LWAVE = .FALSE.
LCHARG = .FALSE.

IBRION = 2
POTIM = 0.1
EDIFFG = -0.02
NSW = 100

3) 准备输入文件(KPOINTS)

1
2
3
4
5
K-POINTS
0
Gamma
1 1 1
0 0 0

4) 准备POTCAR 并提交任务

5) 优化完的结构:

只需要几步优化就完成了。这是因为我们复制的别人的(VASP官网例子)结构,因为这些结构已经被优化过了。即我们优化过程有一个很好的猜测结构。

这给我们的启发就是:在你开始计算的时候,一定要善于利用前人的或者自己前面的计算结构。这会大大节省我们的计算工作量。


末态结构的搭建优化

在上面结构的基础上,我们很容易搭建初末态的结构。仔细观察这个图,我们可以保持H原子的坐标不变,只需要将N原子的坐标修改一下。初始结构里面,三个H原子在xy平面上,

N与H原子z方向的距离为0.3929 Å。如果N原子翻转过去,那么它的坐标将会是:

4.3929 + 0.3929 = 4.7858 Å。

自行将这个图顺时针旋转90°

然后跟前面同样的方法,优化一下就可以了。

过渡态结构的优化

上图中,我们可以发现,过渡态是一个平面的结构,我们可以按照前面末态结构搭建的思想弄一个初始的过渡态结构出来。即直接将N原子的z方向坐标修改成H原子的。如下图:

对比下前面初末态结构的坐标,你会发现Z 方向上被固定住了(最后一列SZ没有选择)

1
2
3
4
+3.8177218605  +3.9717738141  +4.3929644172  T T F
+3.0000000000 +2.5559112549 +4.3926756556 T T F
+2.1822781395 +3.9717738141 +4.3929644172 T T F
+3.0000000000 +3.5000000000 +4.3929644172 F F F

这样做是强制优化过程中,所有的原子都在xy平面上。我们只优化H原面的坐标。原因如下:

这个三角形对应的就是我们从初始结构(IS)搭建过渡态结构(TS)的过程。在初始结构里面,N在三角形最下面,H在右边的点上,当我们人为把N原子放到H平面上的时候,N从下面的点移动到上面,导致N—H 键变短(从1.022 减小到0.944 Å)。但这个变短的过程是我们人为搭建结构所导致的。所以我们要把H原子在xy表面上放开,继续优化一下。参数(ICNAR,POTCAR, KPOINTS)保持不变,提交优化任务,最终我们的过渡态结构如下:

N, H 在同一平面内,N—H 键长为:1.004 Å


频率计算:(过渡态的验证办法之一)

验证一个算完的过渡态是不是真正过渡态的时候,我们需要用到频率计算。因为过渡态只有一个虚频。(自己看书去学习为什么只有一个虚频?)复习下前面我们讲过的频率计算的细节:

1
2
3
4
5
IBRION = 5
POTIM = 0.015
EDIFF = 1E-6
NFREE = 2
NSW = 1

5)此外,我们需要另外加一个无关痛痒的参数:NWRITE = 2 (或者不设置,默认值是2)

大师兄本人喜欢用NWRITE = 0 。但当参数为0,使用Jmol可视化频率计算结果的时候:你会发现所有的原子都是同一个颜色,如下图:

当使用NWRITE的默认值2的时候,则是正常的:这是因为NWRITE = 0 写入的OUTCAR的数据太少,Jmol无法获取元素相关的参数用来可视化。此处,特别感谢群里的:连赞小朋友帮忙分析.

提取频率相关的参数,可以使用下面的命令:

grep cm-1 OUTCAR 提取所有的频率(检查频率结果推荐使用)

grep ‘i=’ OUTCAR 提取所有的虚频 (基本没啥用)

grep ‘f =’ OUTCAR 提取所有的非虚频 (计算零点能推荐使用)

Jmol中所有的虚频都在最下面,大家可以根据前面的grep命令的结果和Jmol的界面对比下。本计算一共有3个虚频: 5.7, 10.2和806.8 cm-1 。前面2个虚频很小,直接不管就可以了。

总结:

算到这里,NH3翻转的过渡态我们就已经讲完了。本节采用的是直接优化的方法获取过渡态的结构,适用于一些非常简单的体系或者你已经非常熟悉的过渡态结构。回到开头我们讲的那句话,过渡态只是一个结构而已,我们需要做的就是如何通过一些计算的技巧和方法来获取这些结构。很多人一提到过渡态就会想到使用VTST编译的VASP程序算NEB(Nudged Elastic Band method)以及CI-NEB (Climbing Image NEB),其实NEB和CI-NEB也仅仅是大家常用的一个方法而已,除此之外,还有很多。

最关键的就是通过结构化学的基础知识来搭建、计算并判断这些过渡态结构背后的物理化学意义,通俗点说就是:算的对不对,你要自己心里有个数。

很多人会说算过渡态要通过频率来验证,其实这一步只是前面一步的副产物而已。只要你的化学基本功底深厚,加以训练,就可以准确判断反应的路径。

本节计算例子的下载: 链接:https://pan.baidu.com/s/1hHGhrOlxMik8NBUqEzolwQ 密码:2dw5

上一节我们知道怎么获取气相分子的熵以及吉布斯自由能的计算。这一节,我们简单介绍一下表面吸附分子熵的计算。

复习基础知识

首先我们应该要知道对一个体系来说,都有哪几部分对熵有贡献。一般熵通过统计热力学计算分子的配分函数来获取:主要有平动,转动,振动,电子以及核这5部分。参考下面的这个链接: http://210.45.168.34:8080/elite/wlhx/jiaocai/C_05.htm 或者找本物化书狂补统计热力学的知识。

1)电子运动的能级间隔很大,除了在几千度以上的高温条件下,电子常常处于基态。一般情况下各激发态对配分函数的贡献都可忽略。

2)原子核的能级间隔极大,在一般的物理和化学过程中,原子核总是处于基态。所以这一项也可以忽略。

3)到现在我们需要考虑的只有平动,转动和振动这三个了。

分子从气相到表面吸附的过程

当分子从气相吸附到表面上,打个不恰当的比喻,就如同图中的苍蝇被粘住在纸板上,因此就不能随便乱飞了,更不用说施展什么空中技巧了。

同样,对分子来说,举个例子,一个非线性的三原子分子,一共有3N个自由度。气相里面为:3个平动,3个转动以及3N-6个振动。但是当它吸附到表面上,由于平动和转动被限制住了,这6个自由度则会转化成振动自由度,也就是在表面上分子有3N个振动模式。如果你可视化表面上吸附分子频率计算的结果,会发现最后6个分别对应的是平动和转动,但它们已经不是气相中的平动和转动了,我们称之为:frustrated translation,frustrated rotation。

(你大妈(平动转动)已经不是十年前的大妈了,你大爷(振动)永远是你大爷)

这两个在表面化学里面是非常令人蛋疼的东西。这个我们后面慢慢讲。至少现在我们知道了,表面的吸附物种有3N个振动方式。所以,平动和转动对熵的贡献也不用管了,直接处理振动即可。(不会算的筒子们小开心下。)

通过振动频率来计算Entropy

振动频率对熵的贡献,随便查阅一本统计热力学的书,都可以找到答案,这里我们的参考书是第十版Atkins的物理化学。

Physical Chemistry 第10版的642页。获取方式:QQ群文件,Books文件下里面就有,自己下载。


图里面圈出来的公式里面,β,h,c,是基本的物理学常数。v (上面加个波浪号)是波数,单位是:m-1。我们就是通过这个公式来计算单个振动模式对熵的贡献。这个很简单,直接把数值带入公式即可,可以摁计算器,也可以用excel等其他工具进行计算。本节,我们分享一个对应这个公式的python小脚本。

通过单个振动频率计算熵的小脚本,图里面的98.485是错的,应该是96.485,kJ/mol 转换为 eV的常数。

脚本解释:

i)VASP得到的频率波数是cm-1,这也是大家所习惯接受的,所以在12行,单位要乘以100先将 cm-1 转化成m-1

ii) 剩下的14-19行是一些基本的物理量,我们默认温度是300K;

iii) 21-27行:是我们定义的一个计算熵的小函数。

x_i 是书里面的 βhcv

pf_l和pf_r 是花括号里面,减号左面和右面的部分

pf 得到的是花括号的结果

然后entropy 是 pf * R

iv)29-31行就是调用函数,然后输出S, TS部分了。这里S的单位是 J K-1mol-1。TS 的单位是我们熟悉的eV。

v)自己对着图片输入一遍就可以得到脚本了,不要问我要现成的。如果自己照着抄的脚本运行出现错误,反复跟图片里面的进行对比,然后进行修改。


熟悉不同的振动频率对熵的贡献

有了这个小脚本,我们就可以计算任何振动频率对熵的贡献了。一般来说,我们VASP的频率计算会有3N个波数,我们需要做的就是挨个算(虚频除外),然后求和即可。这里,大师兄想强调的并不是如何去通过这些频率去算熵(相信到现在你也已经学会了),而是要了解频率大小对熵的贡献程度。下图是我们计算了不同波数(从10到500 cm-1)的振动所对应的熵,以及TS, 单位分别为:J K-1mol-1 和 eV。

通过上图的结果,你会发现:当波数越小的时候,TS越大,随着波数增加,TS越来越小,在490 cm-1, 300K 的时候,连0.01eV都不到。一般来说,正常的振动频率(3n-6,过去以及现在的大爷那部分)对熵的贡献,基本上就可以忽略掉了,因为一般来说正常的频率波数在400-3000cm-1范围之间。

我们现在知道,波数越小,他们对熵的贡献也就越大,而前面我们说的frustrated translation 和 rotation,它们恰恰具有很小的波数,尤其是frustrated translation这部分,大部分都在50cm-1以下,因此他们对熵的贡献不能忽略。然而,它们既有体相声子谱的一些特征,也有气相分子的平动和转动的样子。因此在处理的时候,要小小心心地将它们剥离开来,非常令人蛋疼。正所谓:雄兔脚扑朔,雌兔眼迷离;双兔傍地走,安能辨我是雄雌?

另一方面:由于它们的波数很小,而VASP在小波数这部分的计算很烂,软件计算的结果也会导致很大的误差产生。一般来说大家都直接放弃,忽略这部分了。当然这样的理由很多人不服,明明知道它们波数小,对熵的贡献大,却又不得不忽略。这里,大师兄再给你另外一个解释:摘自Computational Catalysis (RSC) 这本书的第32页。如下:

在这里建议我们将平动单独处理,假设的是分子在二维平面内可以运动。也就是说只损失了一个方向的自由度。经验告诉我们,对不同的过渡金属表面来说,平动部分对熵的贡献可以作为一个常数来处理。所以,又可以不考虑了…. 以上是大家一般都不考虑这部分的贡献的原因。当然,你也可以根据上面的公式计算frustrated translation对熵的贡献,然后其他的按照振动来处理即可。

如果你对这一部分感兴趣,可以参考下这些人的工作:

1) Campbell 这个牛牛的JACS,还有一篇Chemical Review,大家自己找下。

2) 最新的一篇JPCC

3) Tamkin 这个软件的作者以及她发表的一些文章: http://molmod.github.io/tamkin/

4)Norskov 这个牛牛的一本书:Fundamental Concepts in Heterogeneous Catalysis (群文件也有,自己去找)

https://onlinelibrary.wiley.com/doi/book/10.1002/9781118892114

注意: 在这本书里的P33页,不同的吸附位点,构型数目,也会对熵有贡献:如下图:

这个的贡献不是很大,如果你的覆盖度为1/9的时候,R ln((1/9) / (8/9)) = R ln (1/8), 300K的时候,TS = 300 8.314 ln (1/8) / 1000 / 98.485 = 0.053 eV。差不多可以忽略掉。

再提零点能

师兄,既然前面的frustrated的平动和转动我们都忽略了,那么他们的零点能矫正的时候还用不用考虑。

答:小波数范围对零点能的贡献很小,一般来说,考虑或者不考虑它们的贡献,无关紧要。本人一般直接都加上。

思考

1)表面化学中,熵的贡献主要在吸附/脱附的阶段,因为分子损失/获得了大部分的平动和转动,对于表面吸附物种的熵:

i)3N-6的振动部分贡献很小,可以忽略不计;

ii)Frustrated rotation 大约在200-400cm-1左右,也可以和振动部分一起忽略;

iii)Frustrated translation部分,可以近似为二维平动来处理。也可以忽略,一般来说:忽略 是大家普遍的做法。如果理论基础不扎实,只想简单算算,那么你的腰椎间盘还是不要突出为妙。

iv)表面吸附位点的数目对熵的贡献,基本也可以忽略。

2) 而在表面反应的阶段,我们计算反应热和活化能的时候通过下面的公式:

ΔE = E(FS) – E(IS);

Ea = E(TS) – E(IS)

这两个量的计算都是通过相减来得到的,而这个相减的过程中:

i)大部分熵的贡献可以被抵消掉了;

ii)小波数范围内零点能的矫正部分也会被抵消掉,所以零点能校正的时候,大胆地直接地把非虚频的部分直接加起来就可以了。

iii)如果你不信,那么可以通过频率计算一下IS,FS,或者TS的熵,然后计算矫正过的ΔE, Ea,会发现熵的贡献很小很小,所以一般大家都直接考虑零点能,而不考虑熵的贡献,如果有审稿人问你这个问题,用脚本算一下,回答审稿人说影响不大就是了。

3)总结:表面吸附物种的熵,一路忽略,结果就是啥也不算… 不过,不算归不算,但背后的原因或者依据你得搞明白。

前面几节,我们讲解完成了表面上吸附的计算细节。学习到现在,表面吸附的具体流程,如何简化计算,如何判断优化的结构正确与否,这些都是你应该掌握,并且时刻保持思考的事情。尤其是计算存在N种可能性的时候。当你优化出来最稳定的结构,一定要认真琢磨一番,为什么这个结构比其他的构型稳定?是什么因素导致了这个结构最稳定? 从本人多年老司机的经验,简单来说:最稳定的结构从骨子里都透着美的气质,也就是比其他结构看起来更加顺眼。


一旦我们获得了最稳定的结构,就可以进行后面的其他计算了,比如相关的过渡态,与其他分子的吸附能进行对比,与已有的参考文献进行对比,等等。应公众号阅读和分享冠军的要求,我们本节以及后面两节要讲解一下吉布斯自由能的计算,这也是物理化学、化学反应中最常见的一个物理量,尤其是对做电化学的筒子们来说,吉布斯自由能的计算是大家都需要牢牢掌握的。物理化学中,我们学到了:G = H – TS。 而对于G的计算,大家普遍卡在S的这个部分。再次强调,VASP不可以直接算熵,因此我们需要通过频率来自己计算,而OUTCAR里面的那个entropy不是我们物理化学的熵。本节我们简单介绍一下表面吸附物种的频率计算的一些注意事项,以及气相分子的S的获取以及G的计算


1 准备频率计算:

回顾前面我们讲解的乙醇分子的计算例子,然后进行如下的步骤:

1) 修改INCAR的参数

IBRION = 5

POTIM = 0.015 或者0.01-0.02范围的值都可,闲的没事也可以测试下。

NFREE = 2

NSW = 1 1-10000的值均可,没事也可以测试下。

2) POSCAR: 将表面的原子固定住,使用sed命名,或者自己写个小脚本来处理。

3) KPOINTS,一般来说,很多人都使用和优化过程一样的KPOINTS,但从我个人的经验来看,使用更小的K点来计算则是一个更加有效,节约机时的做法。但前提是你要做足够的测试

4) POTCAR保持不变 (废话,等于没说)


2 K点的测试:

回顾前面谁偷了我的机时系列,思考下面的计算,学会分析节约机时,省钱省力省时间。

1) 我们首先使用111, 221, 331, 441进行测试,然后对比振动频率的不同。EDIFFG 对于所有的计算均为 1E-7。

小技巧1): 提取频率的波数,对应的能量,使用 grep cm-1 OUTCAR 这个命令

小技巧2): 只提取虚频,使用 grep f/i OUTCAR

下面是我们测试的结果。H2O为三原子非线性的分子,振动的自由度为:3N - 6 = 3; 也就是有三个振动方式。减去的那6个分别为平动和转动。

对比下不同K点振动波数的区别,都在下面图里了,大家自己看下。

上面图中,前三个为振动的一些能量信息,通过对比可以发现:

i)使用 111 的KPOINTS的时候,我们得到了一个虚频,而其他的KPOINTS没有。

ii)对于振动来说,所有的KPOINTS的结果基本保持不变。


2) 看下不同K值对计算时间的影响。

啥也不想说了,111所花的时间也仅仅是441的一个零头。如果你是个糙哥,111的结果完全够用。如果你是个细心的软男,不敢用111,但又害怕浪费机时,221则是个很好地选择。


3 EDIFFG的测试

这里之所以测试EDIFFG这个参数,是因为它决定了每一个离子步中电子步的步数。收敛条件越苛刻,需要的时间也就越多。

1) 下面是EDIFFG分别为1E-5,1E-6和1E-7的结果, 所有的计算中:KPOINTS为441

通过上面,你可以发现:EDIFFG 从1E-5到1E-7,振动频率变化很小。这个给我们的启示是,1E-5够用了。

2) 我们看一下时间:

K点都为441的时候,1E-5所花的时间是1E-7的一半。


4 思考与测试

i) 通过前面的测试结果, 是不是就可以说我们可以使用gamma点和EDIFFG = 1E-5结合进行计算呢?

2)测试不同的KPOINTS + EDIFFG, 并且与 441 + EDIFFG = 1E-7的结果进行对比,筛选出最佳的节约机时,而又保证结果准确的参数。


5 总结一下表面上频率的计算:

1) POSCAR, 固定表面

2) 修改INCAR: IBRION, POTIM, NFREE, 以及NSW(可有可无)

3) Play with KPOINTS + EDIFFG,节约机时,省时间去做其他的事情,给老板省钱。

4) 如果收敛困难的体系,可以保存优化步骤的WAVECAR用于下一步的频率计算。


6 零点能的计算:

频率计算完成之后,我们便可以计算零点能了。回顾前面乙醇的例子。这里的振动零点能等于前三个振动模式能量之和,(EDIFFG = 1E-7, KPOINTS=441)的结果为:

(462.8 meV + 449.4 meV + 193.1 meV)/ 2000 = 0.55 eV。


7 熵的计算

零点能算完之后,我们还可以继续算以算熵。通过振动频率计算熵的做法我们下一节具体讨论。本节先介绍一下气相分子熵的获取。对于H2O来说,有3个平动和3个转动,而VASP在这方面做的实在是太差了。如果你有兴趣,可以翻阅一下Gaussian对于熵的计算细节:(下载链接http://gaussian.com/thermo/

我们可以根据这个,手动去计算平动和转动对熵的贡献,但相信大部分人不知道怎么去计算。庆幸的是,气相分子的熵值有热力学数据可以查。这里大师兄推荐2个数据库:

1)CRC Handbook

这是本人使用频率最高的数据库了。只要数据库里有的东西,写文章的时候直接引用这本书即可,非常方便。

如果有师弟师妹再问你类似相关的问题,也可以用这本书直接砸砸砸!真可谓是绝世好板砖!


2) NIST 数据库

网址: https://janaf.nist.gov/

如果你看一些原子热力学或者其他相关的计算,你会发现很多人在这个数据库里面获取熵值。很多人不仔细阅读文献,经常给我抱怨说找不到别人怎么算的。是的,你肯定不知道别人怎么算的,因为人家直接查的表,没有算。下面是这个数据库的简单用法:

1)点击红圈中的箭头:

2) 出来的下拉菜单中,点下H,然后找到H2O

3) 点view ,你会得到下面的界面:

-

不同温度下的熵值已经在表中了。直接拿过来用,气相的T*S到手!


8 吉布斯自由能的计算(气相分子)

气相分子的吉布斯自由能为: G = H - TS = U + PV - TS = E_DFT + E_ZPE + nRT - TS (n=1)

i) E_DFT就是我们直接提取的OUTCAR的能量。

ii) E_ZPE 我们计算的气相分子的零点能。(注意:是气相的,不是本节中我们算的表面的零点能)

iii)TS 查表获取

iV) 注意单位要统一。


9 总结:

通过本节,你应该学会

1)怎么计算表面吸附物种的振动;

2)计算零点能;

3)知道频率计算中KPONTS和EDIFFG对计算时间和结果的影响;

4)知道怎么去测试它们的影响,并筛选出合理的参数,节约机时;

5)知道气相分子的熵怎么获取;

6)知道气相分子吉布斯自由能的计算。

7) 关于数据库:如果你是

i) 单身狗:如果有心仪的女生问你怎么算气相分子的熵,请一定不要告诉她这两个数据库。自己查了以后告诉她,没准儿她下次还会来问你。大师兄也只能帮到这儿了,剩下的就看你们自己的发挥啦!

ii)单身的小姐姐或者小妹妹:有了这两个数据库就完全掌握了主动权,剩下的就看个人的心情喽。

iii)老板的话,请直接拿板砖砸!

前面我们搭建了十二个H$_2$O的初始吸附构型。为了获取最稳定的结构,我们需要将这些可能的结构都优化一遍,然后通过吸附能来判断。但是这12个结构,我们都要优化的话,需要很多的机时,尤其是对于那些不太富裕的课题组来说,计算量着实不小。那么我们应该怎么办呢? 回顾:

  • 1)前面几节我们讲的通过简化模型来加快计算的步骤
  • 2)谁偷了我的机时系列。

这里我们可以选择的办法有N个:主要是这两个:1) 减小K点; 2)减小模型。

1) 减小K点的做法

  • i)固定住所有的表面原子(H$_2$O分子除外)
  • ii)使用gamma点进行计算

2) 减小模型的做法

  • i)删去底部的2层原子
  • ii)固定表层的所有原子
  • iii)正常优化进行计算

3) 也可以前面两个方法结合起来进行操作

不管你使用上面什么方法,简化模型这一步的主要目的就是在最快的时间,获取一个最理想的初始构型。当你的初始构型很多的时候,这种办法非常适合作为一个初始的筛选步骤。

师兄,为什么不通过ENCUT来减小计算量?

答:完全可以,只是本人不经常用这种办法而已。个人偏好将INCAR保持不变,以避免ENCUT不同所导致的计算结果错误。


1. H$_2$O 吸附构型的第一步筛选:

这里,大师兄采用的是前面提到的第一个办法,即固定住表面,采用gamma点,然后优化H$_2$O。

  • 1) 准备INCAR,KPOINTS(1x1x1);

  • 2) 批量固定表面:

    1
    sed -i '10,27s/T/F/g'  */*/POSCAR 
  • 3) 将准备好的INCAR和KPOINTS批量复制到各个文件夹中:

    1
    for i in */*; do cp INCAR KPOINTS $i ; done
  • 4) 批量提交任务qsuball.sh


任务结束后,首先查看一下能量信息:

1
for i in */*; do echo $i $(grep '  without' $i/OUTCAR | tail -n 1| awk '{print $7}') ; done

上图,我们发现能量数值相近的差不多,这表明它们很可能是同一个结构。我们把数值sort一下,如下图:

1
for i in */*; do echo $i $(grep '  without' $i/OUTCAR | tail -n 1| awk '{print $7}') ; done | sort -k 2


通过使用p4vasp 批量查看:(p4v //CONTCAR) 你会发现能量为-145.56XXX的都是H$_2$O平行吸附在Cu表面,O在Cu的top位上。


而剩下的,就只有两种结构了:

  • 结构1) type_3/bri 为V型吸附,O在桥式位置上。

  • 结构2)type_3/top为V型结构,O在top位置上。

2. 进一步优化

到现在位置,前面的12个初始结构,可以快速的被筛选成了3个。将这三个结构复制step2的文件夹中。在此基础上,开始正常计算:

  • 1) 将CONTCAR批量复制成POSCAR:

    1
    for i in *; do mv $i/CONTCAR $i/POSCAR ; done
  • 2) KPOINTS变回原来的(4x4x1): sed 批量操作

  • 3) 表面两层原子放开: sed 批量操作

  • 4) 默背一遍提交任务的几个输入文件:INCAR,KPOINTS,POSCAR,POTCAR,脚本,集中注意力思考是否还有没有考虑到的参数或者细节,确保无误后

  • 5) 批量提交任务: qsuball.sh

  • 6)等待结束,查看结果

你会发现两个top的能量也一样了,说明这两个结构很可能优化到一块去了。查看后发现都是平行吸附的结构了。

  • 7) 结构细节:化学上,我们对键长,键角这些信息一定不要放过。测量一下Cu—O 的距离,为: 2.343 $\AA$。

  • 8) 计算吸附能:

与本人4年前的计算结果(-0.17 eV)相差0.01 eV。ACS Catal., 2015, 5, 1027-1036: https://pubs.acs.org/doi/abs/10.1021/cs501698w

本人发表文章支持信息里面的数据。

☆☆☆注意1:

这里吸附能计算的时候,我们取的是最稳定的结构。文献里面提到的也通常是作者所找到的稳定结构。当然,对于不稳定的结构,吸附能也就是顺带的事情了。写在文章里面的话,要标明什么结构对应的吸附能是多少。

☆☆☆注意2:

分子在表面上是物理吸附?还是化学吸附? 通过计算出来的吸附能你应该会判断。如果不知道的话,那么就需要多多翻阅物化书了。


3. 扩展练习:

  • 1) 根据结构和吸附能判断我们的计算结果是不是对的?
  • 2) 查找其他相关文献,将自己的计算结果与文献的进行对比。
  • 3) Step1中12个计算,以及step2中的3个计算,结构和输入文件已经打包,下载链接:https://pan.baidu.com/s/1iyS_nzhI-MuJykabI6rj1w

4. 总结:

通过本节,你应该学会的知识有:

  • 1) 学会如何通过降低工作量来快速筛选不同的初始结构,并最终获取稳定的构型。
  • 2) 巩固和加深对批量操作的理解。
  • 3) 对于更加复杂的结构,其吸附能的基本计算流程要掌握。

结合上一节的内容,下面主要讲解一下H$_2$O分子在Cu(111)表面吸附的初始模型搭建。

1. 计算H$_2$O气相分子的能量

相信大家已经轻车熟路了(老司机)。只将INCAR,KPOINTS,POSCAR列出来,不再详细解释。如下图:

提交任务,结束后,将CONTCAR中的Direct坐标转换为Cartesian,方便后面的表面模型搭建。

怎么转换前面也讲过了,这里也不再啰嗦。


☆☆☆注意:很多人问H$_2$O的分子计算用不用打开自旋的问题?

我的回答是:

  • 1)如果你的体系是开壳层的,那么就打开自旋。
  • 2)如果你的体系是闭壳层的,那么就闭嘴。
  • 3)如果你不知道自己的体系是开壳层的还是闭壳层的,那么就:
    • i) 自己去测试一把,然后对比下,打开与不打开的区别
    • ii)查找文献,看别人怎么做的。

2. 搭建H$_2$O的吸附模型

1) 打开p(3x3)-Cu(111)的slab结构。

先点击左侧的Build,然后点击右上方的text 按钮,进入编辑的部分。如下面的示意图:

注意:上面的只是p4vasp的一个界面,不是p(3x3)-Cu(111)的slab结构。


2) 元素行,原子数目

在元素行,添加O和H,下面一行,添加对应的原子数。

3)添加H$_2$O分子坐标

在坐标的最后面,将前面CONTCAR的Cartesian坐标直接复制过来即可。

4) 结果:

师兄,水分子跑到slab里面去了,这可怎么办? 答:往上平移下即可。

选择O H 以及表面上的任意一个Cu原子。回顾前面的平移操作:

From 那行点 1st 代表的是O原子

To 那行点 last 代表的是表面的4号Cu原子

得到的vector 就是从O到Cu的平移矢量。

☆☆☆注意:

Vector的第三个数字,这个是O和Cu在垂直方向上的距离。我们知道,H$_2$O在Cu上吸附的时候,Cu—O有一定的距离,不妨先给个初始值2A。所以,我们需要把 1.246改成3.246。


小窍门:在下面红色圈出来的部分,可以通过一些快捷的操作来实现原子的选择。

  • i) 在p4vasp 里面可以直接通过写元素符号,选择该元素所有的原子;
  • ii) -10 代表从第一个到第10个原子;
  • iii) all 代表所有的原子;
  • iv) 10- 代表从10到最后的原子;
  • v) 10-25 代表从10到25号原子;
  • vi) 剩下的,自己瞎捉摸去吧,我也不知道了。

上面移动的效果如下:

5) 周期性显示的小技巧

刚开始选择的4号原子在两个格子之间的界面上,本人喜欢将分子放到表面的中间。操作如上图:

  • i) 选择表面上的2个原子: 4 和 5 来定义平移的一个vector。
  • ii) 平移一下O和H原子,效果如下:

当然,在开始的时候,你也可以将H$_2$O分子直接移动到5号原子的表面,大师兄这里做了一些无用的分解操作,主要是为了加深大家的操作印象。


3. 考虑不同的可能性

到目前为止,我们已经顺利搭建了一个H$_2$O分子的吸附模型。把它保存到type_1目录下,名字为POSCAR。

师兄,为什么要建一个type_1的目录?

因为除了这个结构,H$_2$O分子在表面上还有很多其他的可能性。为了保证我们可以获得最稳定的吸附结构,我们就需要将这些可能的结构都算一遍。那么还有其他什么样的可能性呢?

  • 1) 一个H原子指向天空,另外的O—H 键平行于表面。这个结构怎么实现? 回想上一节我们讲到的绕着O—H键旋转的操作。看下图,不解释。

  • 2) 2个H原子都指向天空,只留一个O原子在表面。(剪刀手V)


操作思路: 定义一个穿过O原子的y 轴,然后旋转氢原子。

将这两个结构分别保存到type_2 和 type_3 目录下,名字为POSCAR。按照这样完成的话,你的目录应该是这个样子的:

☆小窍门:p4vasp */POSCAR 打开当前所有目录下的POSCAR,然后:

☆☆☆注意:本操作是在Ubuntu系统下操作的。Windows的用户可能不能实现p4v */POSCAR 这样的操作,下面讲解的搭建模型可能会多花出一些时间。在这里大师兄表示无能为力。


4. 不同的吸附位点

我们前面了解到Cu(111)表面也有很多不同的位点:Top,Fcc,Hcp,和 Bri。所以这三种的初始构型在这四个位点上也要考虑。

Bridge的操作,选择5 和 6, 然后将Vector的数值减小一半,然后平移。

平移完了,将结构保存到bri目录下。

此时不要关闭平移的对话框(这是☆重点1)

点system对应的这个框,选择type_2 的POSCAR(☆重点2)

再点一下move,就完成这个类型的bri结构了。

同样,不要关闭对话框,打开type_3的结构,继续平移,获取对应的bri位结构。

这样操作的话,避免了重复获取平移Vector的那些鼠标点点的操作。点一次,可以继续用于其他结构的操作。希望大家可以掌握这个小窍门,提高自己搭模型的速度。Windows的用户,可以通过将Vector的数值先记下,然后打开另一个结构的时候,直接输入vector就可以了。类似的,我们可以快速搭建完FCC,HCP的结构,最后,你的目录应该是这样子的:

将这些结构复制到服务器里面,准备提交任务。


5. 扩展练习:

  • 1) 完成这12个结构的搭建,彻彻底底地掌握p4vasp平移,旋转操作。
  • 2) 使用MS,VESTA,p4vasp 搭建一个p(4x4)-Pt(111) 的吸附模型。(Bulk,slab的优化过程不用做,只熟悉搭建模型的过程。
  • 3) 思考如何批量准备输入文件,并提交任务?
  • 4) 思考准备的任务输入文件有什么需要注意的地方?

6. 总结:

学完这一节:

  • 1) p4vasp的基本操作,我是没有什么其他可以再教的了。基本上学会这些可以说p4vasp已经入门了,剩下的就是大家没事瞎点点,然后琢磨的事情了。
  • 2)至少含有1-3原子的分子吸附模型你应该会搭建了。原则如下:
    • i)通过基本的化学常识和自己所学的知识判断与表面结合地原子;
    • ii)分析表面不同的吸附位点;
    • iii) 思考可能的吸附方式(躺着、竖着、斜着、歪着等等),并通过平移,旋转来实现。

而最后一点,每思考一种吸附方式,都会使得自己的模型数目激增(这种吸附方式在不同位点上),增加工作量。所以,大家在操作的时候,一定要多看文献,因为前人已经把这蛋疼的路帮你走通了。但这并不意味着别人算的就是对的,大家都是人,也会有疏忽的时候,如果你感觉某个结构应该会更稳定,一定要去尝试一下。

经过前面的学习,相信对于:单原子和双原子分子的吸附,大家已经掌握了怎么搭建模型,以及计算吸附能了。下一节我们要讲的是三原子分子的吸附。在进行三原子分子吸附的计算之前,我想是时候祭出本人使用p4vasp搭结构的一个法宝了。古人语,工欲善其事,必先利其器。所以,本人认为这一节的学习对于以后的结构搭建非常重要,尤其是对于那些使用p4vasp的筒子们来说。说了这么多废话,本节主要内容是手动搭建一个H$_2$O分子,并实现分子的任意旋转。


1. 获取结构

首先我们先知道H$_2$O分子的基本结构。本人的做法是去NIST数据库查找资料。操作如下图:

点击上的链接,或者点击 NIST网址查询得到下面的信息:

在NIST的数据库里面,基本的热力学数据,结构都有了。大家完全可以按照上图中的Cartesian坐标搭建一个H$_2$O的结构,直接复制到POSCAR即可。不过,这个操作我们先缓一下,以后你有的是机会操作。本节主要介绍下如何是从头开始搭建一个H$_2$O分子模型,整个操作过程的学习比搭建模型的结果更加重要。


2. 模型搭建:

1)数据库结构

O—H的键长为0.9578 A,∠HOH = 104.478°。

2)直角形H2O分子

首先,我们在xy平面内搭建一个直角形的H$_2$O分子。O放到原点: 0, 0, 0,H_1沿着x轴, 0.96, 0, 0,H_2 沿着y轴: 0, 0.96, 0

3)修改角度

前面说了,∠HOH = 104.478°, 我们该怎么做呢?

  • i) 可以自己手动算一算,其中一个H原子的坐标,然后更新坐标信息。三角公式大师兄早就忘的一干二净了,暂时跳过。

  • ii)通过p4vasp进行旋转操作。

    对于旋转操作,我们首先要定义一个旋转轴;然后选择旋转的原子,以及角度。在上面H$_2$O的结构里面,很容易就想到,如果∠HOH 从 90°增大到 104.478°, 我们需要以O原子所在的z轴,旋转一个H原子即可,旋转角度为14.478°。


4) p4vasp获取旋转轴的操作:

  • i) 选中O原子;
  • ii) P4vasp界面左上方, Edit—>Rotate Atoms
  • iii) 弹出的界面点击:Get Group
  • iv)Center 那一行,点击1st的按钮;
  • v) Second point那一行,点击1st的按钮;
  • vi) Axis那一行,点击 Z的按钮。

这样我们就定义了一个穿过O原子的z轴。

5) p4vasp选择要旋转的原子和角度:

  • i)把鼠标分别移动要要旋转的H上,通过键盘的space 空格键选择。
  • ii)在上一步Edit –> Rotate Atoms 弹出的对话框中,再点一下: Get Group
  • iii)最后一行,我们选择 -15 作为旋转的角度。

注意1: +15 和 -15 顺时针和逆时针旋转15°。 如果你旋转错了,不要心慌,把数值改下,重新旋转即可。

注意2: 你也可以先选中O和H原子,然后在center 和 second point 那两行都选择1st来获取旋转轴。

点击最下面的: Rotate,效果如下:

由于我们绕着O原子旋转,所以O原子选中或者不选中,旋转操作都对其坐标都没有影响。通过Structure – Measure 测定一下 三个原子的角度。


3. 绕任意轴旋转

上面我们介绍了一下,绕一个单原子的轴进行的旋转操作。而实际的模型调整,搭建的过程中,这种情况并不多,大多时候我们需要绕着一个3D空间里面的轴进行旋转,而不仅仅局限在xyz这样简单的情况。我们知道,两点可以确定一条直线,所以,对于旋转轴来说,我们可以通过两个原子来定义。下面我们讲解一下,H$_2$O分子绕着一个O—H键的旋转操作。

1)定义旋转轴:

  • i) 选择 O 原子 和 H原子;
  • ii) Edit –-> Rotate Atoms;
  • iii) Get group;
  • iv) Center 那一行,点 1st;
  • v) Second point那一行,点 last。

点完之后,下面一行会自动填充我们选择的旋转轴。即从1st原子指向last原子的一个轴。最后一部分,我们写的是60,即绕着O—H键 旋转60°。

选择要旋转的原子,由于前面两个原子已经用来定义旋转轴了,剩下的第三个就是我们旋转操作的对象。选中所有的原子,然后点击 Get group。 下面的图是:H$_2$O 的一个H原子绕着另外两个原子的轴旋转 -60 °的效果。我们绕着一个O—H 键旋转,旋转操作对轴上的原子坐标没有影响,所以大家可以选中这些原子,也可以先通过这两个原子定义旋转轴,进行最后一步操作的时候,取消选择也可以。

2)移动水分子

一般对于分子的计算来说,如果分子在原点附近的话,由于周期性的原因,结构的一部分会进入到另外一个相邻的格子里面,虽然对于计算没有什么影响,但对后面的其他可视化过程(比如,频率计算等)会造成一定的影响。所以本人经常把原子移到格子里面。选择所有的原子,在最后的Vector部分,选择5 5 5, 即在xyz三个方向上都移动5A的距离。点击move,效果如下图,然后保存结构。


4. 扩展练习:

  • 1) 重复本节的所有操作;
  • 2) 复习前面乙醇分子模型的搭建
  • 3) 根据今天所学,随意操作乙醇分子中,原子的旋转,平移等操作,直至熟练位置。

5. 总结:

本节的重点是学会用p4vasp 进行分子的旋转操作,并复习下分子的平移操作。分子的旋转在表面结构模型的搭建过程中,非常重要,熟练掌握这一个技巧,可以极大提高自己搭建模型的合理性和准确性,从而在后面的计算过程中,节约我们的计算时间。如果你有更好的方法,也也可以分享经验给大家。

前面我们讲解了一堆单原子的吸附计算和一些日常的操作,现在是时候来点复杂的操练了。这次,我们拿CO开刷,计算它在Cu(111)表面上的吸附。由于CO比前面的O多了一个原子,复杂性稍微提高了些,但也不是太复杂。作为一个由简单向复杂体系的过渡,是一个很好地例子。

在计算CO的吸附之前,我们首先要了解以下3点:(都是教科书里面的经典内容,下面只列出来,不再细说,如果不懂的话,找本结构化学书好好啃一啃。)

1. CO分子信息:

1)CO的分子结构,如下图:

2) CO的几何结构:(用来搭建吸附和气相的模型)


3) CO与金属的成键方式 (用来搭建吸附模型)

2. 搭建CO吸附的模型:

上一节我们学到了闭着眼操作的一些方式,这一节我们继续闭着眼搭结构。操作流程如下:

1) 将上一节中 C 原子吸附的那几个模型复制过来,进行修改:

  • i) 修改POSCAR中的元素(第6行)和原子数目部分(第七行),添加O原子。

  • iii)在POSCAR的尾部,添加O原子的坐标:

可以看到上图中,最后两行的坐标一模一样,这是因为CO在表面上吸附的时候,我们假定的是直立吸附。所以C 和 O原子的坐标在x和y方向上一样。不同的区别在z方向上,即CO的键长。CO如果吸附在表面上,肯定会和表面原子有作用,也就是所谓的活化,如果原子被活化了,那么CO键就会被削弱,具体体现在键长上。与气相的键长相比,表面上CO的键长数值更大一些。前面我们查数据库得到CO键长为:1.138A.这里我们不妨设置成1.2A。根据这些,我们设置O原子的坐标,如下:

同样,我们可以对其他吸附位点的坐标进行类似的修改,结果如下图:

3. 背一遍VASP的输入文件,检查还有那些需要修改的:

  • i) INCAR: 跟之前保持一致;
  • ii)KPOINTS:跟之前保持一致;
  • iii) POSCAR: 已经修改完毕;
  • iV) POTCAR:提交任务的脚本里面自动生成
  • V)提交任务命令: qsuball (前面已经讲解过了)

4. 提交任务:

上图,我们删除了slab的计算,因为前面我们已经计算过了,没有必要浪费机时再算一遍。

5. 思考下吸附能的计算公式:

前面我们还忘记了CO的气相结构优化。现在我们回忆一下前面所讲的气相分子的优化:

1) 气相分子的结构模型搭建。

直接将一个Cu(111)表面上吸附的POSCAR拿过来修改一下即可:

2) INCAR:直接拿Cu(111)吸附的修改下。

对于气相分子的优化来说:

1
2
ISMEAR = 0 
SIGMA = 0.01

3) KPOINTS: gamma 点即可。

4) 这里我们没有批量提交,手动运行一下:pospot.sh脚本,生成对应的POTCAR。

5) 使用提交单个任务的脚本: qsub 提交任务。具体操作如下图:

6. 扩展练习:

  • 1)完成CO吸附的计算;
  • 2)进一步熟练简单模型的闭着眼操作方式;
  • 3)复习分子气相结构的优化过程。

7. 总结:

本节,通过CO的吸附模型搭建,带领大家走出简单的单原子吸附,开始逐渐接触复杂的吸附计算。如果前面内容掌握了,本节就是一个水到渠成的事情。在多原子分子的吸附计算中,首先我们要知道分子的电子和几何结构,分子哪一部分(这里的C原子)和表面成键。在本节的例子中,如果你不知道C和金属作用,你还需要计算Metal-O-C这样的结构。可能还会计算CO横着吸附的结构,任务无形中就会增加一倍或者更多,从而造成机时的浪费。退一步来说,如果你真的不知道吸附是以什么方式进行的,想尝试N种初始的结构,我的建议是:把slab的原子全部固定住,然后用gamma点算一下它们的吸附能,先大体上判断一下,把那些吸附特别强的结构筛选出来,用作下一步的计算。