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$):

$$ E_a = E(\textrm{TS}) – E(\textrm{IS}) \
\Delta E = E(\textrm{FS}) – E(\textrm{IS}) $$

考虑零点能校正:

$$
\begin{align}
E{a(\textrm{ZPE})} &= E_{\textrm{ZPE}}(\textrm{TS}) – E_{\textrm{ZPE}}(\textrm{IS}) \
&= E(\textrm{TS}) + \textrm{ZPE(TS}) – E(\textrm{IS}) – \textrm{ZPE(IS}) \
&= E_a + \textrm{ZPE(TS)} – \textrm{ZPE(IS)}
\end{align}
$$

同理:

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

两个处理方式:

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

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

效果是一样的

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

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


例子A:

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

$$
\ce{CH_3CH_2OH + Cu(111) -> Cu(111)-CH_3CH_2OH}
$$
此时,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 原子进行振动。

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


还有一点要注意的是:ISTSFS中,所固定和放开的原子必须一致!!!

不要在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 虚频

(这个命令只展示一次,后面测试结果全部用下面的表格)

EDIFFNthWNUnitEnergyUnit
42510.183435 cm-11.262585meV
2643.891375 cm-15.441839meV
2761.646745 cm-17.643225meV
5256.773094 cm-10.839757meV
2628.045092 cm-13.477149meV
2752.28955 cm-16.48308meV
6256.132037 cm-10.760276meV
2618.62878 cm-12.309675meV
2765.657693 cm-18.140519meV
7256.03848 cm-10.748676meV
2619.115009 cm-12.36996meV
2765.726299 cm-18.149025meV
8256.036457 cm-10.748425meV
2619.156334 cm-12.375084meV
2765.711053 cm-18.147135meV

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 虚频:

ENCUTNthWNUnitenergyunit
400256.132038 cm-10.760276meV
2618.62878 cm-12.309675meV
2765.65769 cm-18.140519meV
500240.489004 cm-10.060629meV
2516.7313 cm-12.074418meV
2623.491849 cm-12.912619meV
2766.062971 cm-18.190767meV
600240.288462 cm-10.035765meV
254.853849 cm-10.601801meV
266.324527 cm-10.784142meV
2752.075614 cm-16.456556meV
700250.747158 cm-10.092636meV
266.175845 cm-10.765708meV
2732.540704 cm-14.034535meV
800250.289332 cm-10.035873meV
267.356367 cm-10.912074meV
2726.44812 cm-13.27915meV

增加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, AccurateHigh

3.3.1 零点能变化甚微

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

3.3.2 虚频

PRECNthWavenumberUnitEnergyUnit
N256.132037 cm-10.760276meV
2618.62878 cm-12.309675meV
2765.657693 cm-18.140519meV
A255.493108 cm-10.681059meV
2610.478293 cm-11.299143meV
2765.861524 cm-18.165791meV
H240.552665 cm-10.068522meV
2516.954329 cm-12.10207meV
2629.95621 cm-13.714098meV
2768.287802 cm-18.466611meV

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虚频

POTIMNthWavenumberUnitEnergyUnit
0.0012521.746998 cm-12.696285meV
2639.656607 cm-14.916794meV
2756.54809 cm-17.011072meV
0.005240.675194 cm-10.083713meV
2512.785399 cm-11.585188meV
2636.264179 cm-14.496187meV
2761.287932 cm-17.598738meV
0.015255.626635 cm-10.697614meV
2623.117605 cm-12.866219meV
2759.01173 cm-17.316525meV
0.020256.132037 cm-10.760276meV
2618.62878 cm-12.309675meV
2765.657693 cm-18.140519meV
0.050240.327964 cm-10.040662meV
2523.575077 cm-12.922938meV
2624.337978 cm-13.017526meV
2795.239085 cm-111.80815meV
0.100231.574379 cm-10.195198meV
245.865037 cm-10.727172meV
2520.89784 cm-12.591003meV
2645.241296 cm-15.609208meV
27168.238904 cm-120.85897meV

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 虚频:

KPOINTSNthWavenumberUnitEnergyUnit
1256.132037 cm-10.76028meV
2618.62878 cm-12.30968meV
2765.657693 cm-18.14052meV
2255.726298 cm-10.70997meV
2619.193163 cm-12.37965meV
2764.645347 cm-18.015meV
3254.826986 cm-10.59847meV
2619.138554 cm-12.37288meV
2764.83042 cm-18.03795meV

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}$ 或者水分子的计算作为例子。那样的简单例子不具有很强的代表性,当你把一个相对复杂的东西搞明白了,这些小分子的计算,就如同小儿科一般了。