现在位置: 首页 / 未分类 / 正文

Ex27 乙醇分子的振动频率计算(五)

通过前面读者的反馈，本节前面主要讲述一下零点能相关的注意事项以及频率计算的时间成本相关的知识。

1 参考书籍：

以下所提到的书籍都可以在QQ群文件和百度网盘里面下载：http://pan.baidu.com/s/1bKHMjG

1.1 Density Functional Theory：A Practical Introduction

大师兄极力推荐的菜鸟入门书！！！

1）不要网上随便下载（很多都是野鸡版，公式符号不全）

2）不要看中文版的；去看原汁原味的，困了就歇会，慢慢坚持看；

3）大师兄QQ群下载或者点击百度网盘链接；

4）第5章专门讨论的频率相关的问题！！！

1.2 Atkins的物理化学第9版

前面几节我们都是用的这一部分作为参考，忘了的童鞋们返回继续重修；这本书网上已经可以下载第十版的pdf文件了。

1.3 Jens. K. Norskov 的 Fundamental Concepts in Heterogeneous Catalysis

2 频率计算的细节问题

当你阅读完前面推荐书中的频率部分后，应该会明白很多。本节的第二部分我们主要讨论下频率计算的一些大家有顾虑的细节问题：

2.1 怎么进行零点能校正？

以乙醇的计算为例：

A）首先结构优化完毕后我们会得到分子的能量：$E_{0}$ （OSZICAR中的E0）

B）频率计算后我们会得到分子的零点能：$ZPE$

C）零点能校正之后的乙醇分子能量为： $E_{ZPE} = E_0 + ZPE$

D）A 和 B 得出的结果直接相加即可，不要想太多。

2.2 怎么计算过渡态和反应热的零点能校正？（菜鸟暂且跳过）

对一个反应：IS —> TS —> FS

IS: Initial State 反应物

TS: Transition State 过渡态

FS: Final State 产物

记住这几个定义，本书后面 IS, TS, FS 用的时候多着呢！

1）优化反应物 IS 和产物 FS 的结构，获得能量：$E(\textrm{IS})$, $E(\textrm{FS})$;

2）对反应物和产物进行频率计算，获得各自的零点能：$\textrm{ZPE(IS)}, \textrm{ZPE(FS)}$。

3）搜索过渡态，获得结构和能量 $E(\textrm{TS})$;

4）过渡态频率分析，获得零点能 $\textrm{ZPE(TS)}$

不考虑零点能的反应能垒 ($E_a$) 和反应热 ($\Delta E$)：

考虑零点能校正：

同理:

$ \Delta E(\textrm{ZPE}) = \Delta E + \textrm{ZPE(FS)} – \textrm{ZPE(IS)} $

两个处理方式：

a. 先获取未校正的结果，然后把零点能各自相减；

b. 先将各个物种进行零点能校正，然后在计算反应能垒或者反应热

效果是一样的

2.2 频率计算的时候，是不是体系中所有的原子都放开？

不一定；这取决于你的体系，以及你想要关注的部分:

例子A:

乙醇在Cu(111)表面上的吸附，计算吸附热的零点能校正

此时，Cu(111)表面我们在计算频率的时候是要固定住的！只振动乙醇分子即可。

例子B:

计算CO在Cu(111)表面上的吸附：

同A，固定Cu(111) 表面，如果你只关心 CO 在垂直表面上的振动，那么CO的 $x, y$ 方向便可以固定住，在坐标后面为：F F T

例子C:

苯酚在Cu(111)表面上 O—H 键断裂活化能的零点能校正：

在这里我们拿苯酚作为例子，很多时候，计算的对象比较大，全部频率优化非常耗时，那我们就得选择性地固定住一部分，只关心关键的局域部分。这个例子中我们主要讨论零点能对 O—H 键断裂活化能的影响，因此我们可以把苯环的部分固定住，只放开 O 和 H 原子进行振动。

当然啦，这是一个简化的计算，肯定不如全部振动的结果好，但迫于计算时间的压力，这是一个折中的好办法。此外，由于在反应中苯环部分变化不大，所以即使我们全部放开，它对反应热和反应能垒的影响也很小，零点能部分相减的时候抵消了。在这里你要深刻体会到，意识到，把握到相对这两个字在理论计算中的巨大作用！！！

还有一点要注意的是：IS，TS和FS中，所固定和放开的原子必须一致！！！

不要在IS中固定苯环，而在TS或者FS中放开， 这样计算出来的结果只能用脑残来形容了。(ps. 本人犯的这种错误太多啦，已经不能用脑残来形容了！)

3 影响频率计算的因素测试:

前面提到过：ENCUT, KPOINTS, POTIM, EDIFF都会对频率计算产生影响。频率计算很多时候会有小的虚频出现，通俗点说就是数值噪音。我们可能会想到，凡是可以提高精度的办法是不是就可以降低噪音了呢？写到这里，大师兄突然想起初中的一个老师说过一句话：好记性不如烂笔头，意思是让我们多写，这样记得会比较快。还有一句名言：纸上学来终觉浅，绝知此事要躬行！学习 VASP，亲自上手实践及其重要，看遍官网，翻烂本书不去行动也是白搭，好听点说只能是纸上谈兵。

下面的测试主要从零点能和虚频，和时间三个方面来判断：测试的命令为：

获取虚频命令：

1
 
grep 'f/i'  */OUTCAR | awk '{print $1 "\t " $2 "\t" $8 "\t " $9 "\t" $10 "\t" $11}'

获取时间命令：

1
 
grep Elapsed */OUTCAR | sort -n

获取零点能：

1
 
for i in * ; do echo $i $(cd $i ; fsum ; cd $OLDPWD);done  | sort –n

3.1 EDIFFG测试

3.1.1测试

$4, 5, 6, 7, 8$ 分别代表 EDIFF=1E-4 到 -8 的测试任务；

1
 
for i in * ; do echo $i $(cd $i ; fsum ; cd $OLDPWD); done

零点能的计算结果：

EDIFF ZPE Unit
4 2.14344 eV
5 2.12494 eV
6 2.11698 eV
7 2.11696 eV
8 2.11696 eV

a. 可以看出来在1E=-6以后，零点能变化基本为 0了；我们可以用这个设置；

b. 虽然你测出来一个很好的参数值，很不幸，这是默认的。

c. 从时间上考虑，本人会选择 EDIFF= 1E-5，这比 1E=-6 快了将近一倍；

d. $5$ 和 $6$ 之间的误差为：$ 0.008/2.117 = 0.38\% $；

e. 此外，单个的零点能貌似没有什么用处，一般都是两个结构的零点能相减。这个操作也会抵消一部分的计算误差。

3.1.2 虚频

（这个命令只展示一次，后面测试结果全部用下面的表格）

EDIFF Nth WN Unit Energy Unit
4 25 10.183435 cm-1 1.262585 meV
26 43.891375 cm-1 5.441839 meV
27 61.646745 cm-1 7.643225 meV
5 25 6.773094 cm-1 0.839757 meV
26 28.045092 cm-1 3.477149 meV
27 52.28955 cm-1 6.48308 meV
6 25 6.132037 cm-1 0.760276 meV
26 18.62878 cm-1 2.309675 meV
27 65.657693 cm-1 8.140519 meV
7 25 6.03848 cm-1 0.748676 meV
26 19.115009 cm-1 2.36996 meV
27 65.726299 cm-1 8.149025 meV
8 25 6.036457 cm-1 0.748425 meV
26 19.156334 cm-1 2.375084 meV
27 65.711053 cm-1 8.147135 meV

a. 改变收敛标准，并没有消除虚频，测试 $5$ 中，稍微降低了些，其他情况下，最大的仍然都再 1500 px-1 左右

b. $6$, $7$, $8$ 中的虚频基本一致，说明增加收敛标准，并没起到什么好的效果。

3.1.3 时间

EDIFF Time Unit
4 1702.463 s
5 2295.233 s
6 4298.571 s
7 5481.423 s
8 5998.518 s

复习下谁偷走了我的机时中EDIFF的因素；增加EDIFF，计算时间变长了。（离子步数一致，每个离子步的收敛步数增多了） 一般来说EDIFF采用默认值或者1E-5都可以满足情况，任务比较繁重的时候议用1E-5。本节在后面的测试中，都采用 VASP 的默认EDIFF值进行（1E-6）

3.2 ENCUT测试

3.2.1 零点能

ENCU ZPE Unit
400 2.11698 eV
500 2.11448 eV
600 2.11880 eV
700 2.11971 eV
800 2.11928 eV

从表中可以得出结论：增加ENCUT,对零点能的影响很小；

在计算的过程中，ENCUT要保持一致，一但选定后就不要再轻易改动，这个表的数据告诉我们，即使改动了ENCUT，对零点能的影响也很小，所以ENCUT的影响可以忽略。

3.2.2 虚频：

ENCUT Nth WN Unit energy unit
400 25 6.132038 cm-1 0.760276 meV
26 18.62878 cm-1 2.309675 meV
27 65.65769 cm-1 8.140519 meV
500 24 0.489004 cm-1 0.060629 meV
25 16.7313 cm-1 2.074418 meV
26 23.491849 cm-1 2.912619 meV
27 66.062971 cm-1 8.190767 meV
600 24 0.288462 cm-1 0.035765 meV
25 4.853849 cm-1 0.601801 meV
26 6.324527 cm-1 0.784142 meV
27 52.075614 cm-1 6.456556 meV
700 25 0.747158 cm-1 0.092636 meV
26 6.175845 cm-1 0.765708 meV
27 32.540704 cm-1 4.034535 meV
800 25 0.289332 cm-1 0.035873 meV
26 7.356367 cm-1 0.912074 meV
27 26.44812 cm-1 3.27915 meV

增加ENCUT可以减小虚频；从 $400$ 增加到 $800$ 后，最大的虚频波数从 $65~cm^{-1}$ 减小到 $26~cm^{-1}$， 但 A 中的结果是不是要求我们算频率的时候，必须要把 ENCUT增大呢？ 答：没必要！！

因为前面零点能变化甚微，且虚频从 $65~cm^{-1}$ 减小到 $26~cm^{-1}$, 只能说是在误差范围内变化，我们是可以承受的，此时ENCUT在减小虚频中的作用，大家记住这一点就可以了。此外，增加ENCUT也会使得计算时间变长（见下图）。但是如果虚频大于$100~cm^{-1}$，且不是反应路径的时候，你就得小心了。

3.2.3 时间

命令：

1
 
grep Elapsed */OUTCAR | sort –n
ENCUT Time Unit
400 1832.405 s
500 2389.045 s
600 3019.433 s
700 4285.491 s
800 11943.02 s

3.3 PREC

这里我们设置了三个参数：Normal, Accurate和High

3.3.1 零点能变化甚微

PREC ZPE Unit
N 2.11698 eV
A 2.1147 eV
H 2.11363 eV

3.3.2 虚频

PREC Nth Wavenumber Unit Energy Unit
N 25 6.132037 cm-1 0.760276 meV
26 18.62878 cm-1 2.309675 meV
27 65.657693 cm-1 8.140519 meV
A 25 5.493108 cm-1 0.681059 meV
26 10.478293 cm-1 1.299143 meV
27 65.861524 cm-1 8.165791 meV
H 24 0.552665 cm-1 0.068522 meV
25 16.954329 cm-1 2.10207 meV
26 29.95621 cm-1 3.714098 meV
27 68.287802 cm-1 8.466611 meV

PREC同样对消除虚频不管用！！！！

3.3.3 时间

Prec Time Unit
N 4298.571 s
A 7596.994 s
H 7962.539 s

通过前面的虚频和零点能，结合此时的时间，我们可以得出结论：

PREC = Normal 完全够用了。

3.4 POTIM 参数

3.4.1 零点能：

POTIM ZPE eV
0.001 2.17125 eV
0.005 2.12616 eV
0.015 2.11769 eV
0.020 2.11698 eV
0.050 2.11671 eV
0.100 2.12481 eV

POTIM太小的时候，对零点能影响很大。（$0.001$ 和 $0.005$，$0.015$ 对比）

3.4.2虚频

POTIM Nth Wavenumber Unit Energy Unit
0.001 25 21.746998 cm-1 2.696285 meV
26 39.656607 cm-1 4.916794 meV
27 56.54809 cm-1 7.011072 meV
0.005 24 0.675194 cm-1 0.083713 meV
25 12.785399 cm-1 1.585188 meV
26 36.264179 cm-1 4.496187 meV
27 61.287932 cm-1 7.598738 meV
0.015 25 5.626635 cm-1 0.697614 meV
26 23.117605 cm-1 2.866219 meV
27 59.01173 cm-1 7.316525 meV
0.020 25 6.132037 cm-1 0.760276 meV
26 18.62878 cm-1 2.309675 meV
27 65.657693 cm-1 8.140519 meV
0.050 24 0.327964 cm-1 0.040662 meV
25 23.575077 cm-1 2.922938 meV
26 24.337978 cm-1 3.017526 meV
27 95.239085 cm-1 11.80815 meV
0.100 23 1.574379 cm-1 0.195198 meV
24 5.865037 cm-1 0.727172 meV
25 20.89784 cm-1 2.591003 meV
26 45.241296 cm-1 5.609208 meV
27 168.238904 cm-1 20.85897 meV

POTIM太大的时候，会搞出来超大号的虚频（$0.050$ 和 $0.100$）！

3.4.3 时间

POTIM Time Unit
1 1530.591 s
5 2065.911 s
15 3845.977 s
20 4298.571 s
50 4829.897 s
100 5191.23 s

POTIM 越小，单个离子步收敛的越快，也就是需要更少的电子步！

综合前面的三个因素：POTIM太小或者太大都不好， POTIM = 0.015或者0.020 是很好的选择。 0.015是默认值。

3.5 POINTS

这里用 $1$, $2$, $3$ 来分别代表 $1 \times 1 \times 1$， $2 \times 2 \times 2$ 和 $3 \times 3 \times 3$ 的 K 点。

3.5.1 零点能：

KPOINTS ZPE Unit
1 2.11698 eV
2 2.11678 eV
3 2.11689 eV

3.5.2 虚频：

KPOINTS Nth Wavenumber Unit Energy Unit
1 25 6.132037 cm-1 0.76028 meV
26 18.62878 cm-1 2.30968 meV
27 65.657693 cm-1 8.14052 meV
2 25 5.726298 cm-1 0.70997 meV
26 19.193163 cm-1 2.37965 meV
27 64.645347 cm-1 8.015 meV
3 25 4.826986 cm-1 0.59847 meV
26 19.138554 cm-1 2.37288 meV
27 64.83042 cm-1 8.03795 meV

3.5.3 时间

KPOINTS Time Unit
1 4298.571 s
2 8997.709 s
3 13480.847 s

综合前面三点： Gamma点足矣！

4 扩展阅读：

4.1 下载本节的练习，按照本节的顺序操作分析；

4.2 阅读 VASP 官网中关于原子和分子的例子，去尝试回答里面的问题；

https://cms.mpi.univie.ac.at/wiki/index.php/Atoms_and_Molecules

5. 总结：

分析完前面的测试，我们总结一下高效频率计算的关键点：

5.1 IBRION = 5 （告诉 VASP 我们要算频率）

5.2 POTIM = 0.015

5.3 NFREE = 2

5.4 ENCUT和原来一样

5.5 PREC = Normal

5.6 EDIFF = 1E-5 或者 1E-6

5.7 KPOINTS Gamma点即可。

本节内容太多了，原因在于，这是菜鸟篇的最后一节！！！

大师兄本人亲自做了很多测试，原因只有一个，告诉大家如何把握计算结果与时间的关系。在最短的时间获取最多的有价值的结果，这是计算化学的一个精髓所在。体现了你对时间和生命的尊重，也体现了你对高效率生活的追求。在菜鸟的时候，可以尽情地去测试。

大师兄之所以用乙醇这个大分子去给大家展示，很多人直接用 $\rm{O_2}$，$\rm{CO}$ 或者水分子的计算作为例子。那样的简单例子不具有很强的代表性，当你把一个相对复杂的东西搞明白了，这些小分子的计算，就如同小儿科一般了。