本文章为原创,版权归作者刘锦程所有,文章转载请先取得作者的同意,非常欢迎转发文章链接!严禁以任何方式挪用本文内容,用于以盈利为目的各种活动。部分内容参考一个人就是一个叠加态的博客,
http://blog.sina.com.cn/s/blog_b364ab230102vghk.html,
原理
过渡态(transition state)搜索,优化到一阶鞍点。是MEP路径上能量最高的点(一阶鞍点)。在过渡态位置上,仅在沿着反应路径上是能量极大点,而在于之正交的其他所有方向上都是能量极小点。
NEB是chain-of-states搜索方法的一种,经过长时间的发展形成现在的准确性好,稳定性好的CINEB方法。做法就是在反应物到产物之间插入一系列结构,共插入P-1个,反应物编号为0,产编号物为P。不同的是优化不是对每个点孤立地优化,而是优化一个函数,每一步所有点一起运动。
这P+1点保持第一和最后一个点不动,其他的点受到两个力的作用:1. 来自势能面的力;2. chain方向上的弹簧力。
第一项,相当于让一串珠子上的每个点都在势能面上向着最稳定的方向上移动。
第二项,含有弹簧力常数,相当于每个珠子之间的拉力和排斥力,避免珠子之间的距离过短或者过长。
传统PEB方法(plain elastic band)方法遇到的问题:
问题一:反应路径偏离minimum energy path(MEP),在势能面拐点的地方会出现两个珠子距离过远的问题,如下图,左图。
问题二:难以找到最高点,需要插入非常多的点,浪费计算资源。
NEB方法,通过nudge过程,解决了问题一:
•弹簧力垂直于路径的分量就被投影掉,
•弹簧力平行于路径的分量完全保留;
•势能力在路径方向上的分量被投影掉,
•势能力垂直于路径上的分量将会引导结构点地正确移动。
这样(1) 每个点在平行于路径切线上的受力只等于弹簧力在这个方向分量,(2) 每个点在垂直于路径切线方向的受力只等于势能力在此方向上分量。这样优化收敛后结构点就能正确描述真实的MEP(Minimum
energy pathway),矛盾得到解决。
以上过程讲述的是VASP原版中使用的NEB。但是此方法仍有很多致命的缺点:
(1)要插入足够多的点才可能找到近似的过渡态,计算资源消耗极大。
(2)因为珠子不可能到鞍点的顶端,计算出的过渡态能量总是被低估的,又因为本来LDA,GGA泛函就经常低估过渡态能量,所以误差叠加,此方法不太准确。
所以后来,G. Henkelman(VTST的开发者)又发展出了CI-NEB方法,它专门考虑到了定位过渡态问题。CI-NEB与NEB的关键区别是能量最高的点受力的定义,在CI-NEB中这个点不会受到相邻点的弹簧力,避免位置被拉离过渡态,而且将此点平行于路径方向的势能力分量的符号反转,促使此点沿着路径往能量升高的方向上爬到过渡态。http://sobereva.com/44
CI-NEB方法只需要很少的点,比如包含初、末态总共5个甚至3个点就能准确定位过渡态,是最有效率的寻找过渡态的方法之一。注意: 在网上看到一些说法CI-NEB插点越多越好,这是错误的结论!!由于能量最高的点可以自动爬坡,所以有的时候插一个点也可以精确的找到过渡态的位置,大多数时候插3到4个点已经完全能应付正常的过渡态计算需要。(具体插入多少点合适需要用dist.pl脚本判断)
CI-NEB过渡态计算步骤,练习一:O在Au(111)上的迁移
我们用VTST中CI-NEB做示范,(首先要确保编译VASP的时候安装了VTST,否则用的是VASP自带的NEB,非常糟糕)。
步骤一:优化初态和末态结构。FCC位和HCP位。
步骤二:用dist.pl(VTST脚本)检查两个优化后结构(两个CONTCAR)的相似程度(每个对应原子的初末态距离的平方和,再开根号)。
1 | dist.pl ini/CONTCAR fin/CONTCAR |
若返回值小于5埃,一般可以进行下一步。我的情况:
1.73920306203566
注意:这里如果数值非常大,要检查初态和末态的原子顺序是不是一一对应的!
步骤三:nebmake.pl 插点
插点的数目取决于前面dist.pl的返回值,一般插点数目可取(dist.pl返回值/0.8)。
1 | nebmake.pl ../is/CONTCAR ../fs/CONTCAR 1 |
这个例子插一个点就够了。(好多人被网上的有毒的教程误导,认为插一个点不对,这其实完全没有问题,这里只插一个点和插三个点的结果完全一样。)
返回提示:
1 | filetype1: vasp5 |
代表插点成功,并自动生成了./00 ./01 ./02三个文件夹。
00表示初态,里面放的是../is/CONTCAR, 02表示末态,放的是../fs/CONTCAR, 01是插入的点。三个文件夹里面的文件名称都是POSCAR。
步骤四:初末态对应的OUTCAR复制到对应的文件夹中,以便后续数据分析。
1 | cp ../is/OUTCAR ./00/ |
步骤五:检查插入点的合理性。(可跳过此步)
输入命令
1 | nebmovie.pl 0 |
参数0表示用POSCAR生成xyz文件;还可取1,为用CONTCAR生成。返回值:
1 | Using POSCARs to generate movie |
movie.xyz 直接拖入Jmol可以看所有结构。
1 | 注:用nebmovie.pl,没有直接生成movie.xyz是因为从官方主页下载的脚本默认当使用nebmovie.pl后面不带参数或者参数为0的时候(即使用POSCAR产生xyz),不输出movie.xyz。 你想让脚本自动生成movie.xyz而不是自己去合并各个文件夹里面的xyz文件,需要自行修改nebmovie.pl文件。很简单,把倒数第二个if语句整个用#注释掉或者直接删掉即可。 |
步骤六:准备INCAR
过渡态也是结构优化的一种,所以可以在结构优化的INCAR基础上修改:
1 | (1)EDIFF=1E-7,这个参数非常关键!过渡态对力计算精度要求极高,更精准的电子步收敛会有更精准的力,可以加速收敛。好多过渡态计算不收敛原因都是因为力的精度不够。(注:粗收敛可以适当放宽到1E-5) |
VTST包含的优化算法:
常用的是
IOPT = 1, 2 (适合精收敛)
IOPT = 7 (适合粗收敛)
步骤七:准备POTCAR和KPOINTS
同结构优化直接复制就行
步骤八:提交任务
注意用的核的个数必须整除插点的数量,(为保证最大的并行效率节点数最好等于差点的个数)。
步骤九:跟踪检查收敛结果
1 | nebmovie.pl 1 |
自动生成一个movie.xyz文件,导入到Jmol可以看反应路径上所有点的结构
计算过程中,可以用nebef.pl
或者nebefs.pl
命令(不带参数)查看计算收敛情况:
后三列: [最大原子受力]、 [能量]、 [相对初态的能量]
nebef.pl
: (有多少个点就会显示几行,这是插了4个点的算例)
当所有插点的最大原子受力都<|EDIFFG|时,计算收敛。若中间插了多个点,则所有点都要<|EDIFFG|才能收敛。
步骤十:频率计算验证虚频(下节课)
步骤十一:nebresults.pl总结结果
nebresults.pl也是VTST脚本,用来总结neb结果
输出如下:
1 | Unziping the OUTCARs ... done |
nebresult.pl做的事情如其所输出表明的,完成了nebbarrier.pl,
nebspline.pl, nebef.pl, nebmovie.pl, nebjmovie.pl, nebconverge.pl还有对各文件夹中的OUTCAR打包压缩。生成了很多文件。其中mep.eps是以dist.pl距离为横坐标,能量为纵坐标画出的能势垒图,如下:
注:
nebresult.pl会把所有OUTCAR都打包成.gz的文件,如果不想打包可以用如下命令解压。
1 | gunzip 0*/OUTCAR.gz |
生成的vaspgr文件夹内是各个插点结构的收敛图(同样可用EPS/PS
viewer打开),如下(我的只插入一个结构,所以只有一个图):
转载请注明来源,欢迎对文章中的引用来源进行考证,欢迎指出任何有错误或不够清晰的表达。