||
VASP可以加入VDW修正,但是里面的情况非常多,对于不同的体系该如何选择不同的修正呢?把这个问题先拿出来慢慢的补充讨论吧。
首先先引用侯老师在小木虫的一段话[1]:
有两大类的范德华斯作用修正:
1)基于半经验的,包括了D2, D3, D3-BJ, TS, TS+SCS等等,这些都是在常用的交换关联泛函比如PBE的基础上,考虑色散力的作用,在PBE计算出来的总能基础上增加了额外半经验项,这一项需要设置一些参数,至于是哪种半经验公式,由IVDW的设置来决定,具体的参考VASP的手册。
https://cms.mpi.univie.ac.at/vas ... ection_methods.html
2)另一类就是vdW-DF,也就是修改交换关联项来考虑范德华修正。包括你所提到的optPBE-vdW, optB88-vdW, and optB86b-vdW等。它们的设置由修改GGA赋值,以及额外相关的一些参数来设置。 这些在vasp手册上有明确的说明。
https://cms.mpi.univie.ac.at/vas ... undqvist_et_al.html
那我们就基于这两个大的方面来讨论吧。
第一步还是把vasp官方的说明给出来。
IVDW IVDW = 0 | 1 | 10 | 11 | 12 | 2 | 20 | 21 | 202 | 4 Default: IVDW = 0 Description: This tag controls whether vdW corrections are calculated or not. If they are calculated IVDW controls how they are calculated. Popular local and semilocal density functionals are unable to describe correctly van der Waals interactions resulting from dynamical correlations between fluctuating charge distributions. A pragmatic method to work around this problem is to add a correction to the conventional Kohn-Sham DFT energy : The correction term is computed using some of the available approximate methods. The choice of vdW method is controlled via the following tags: IVDW = 0 no correction IVDW = 1|10 DFT-D2 method of Grimme (available as of VASP.5.2.11) IVDW = 11 zero damping DFT-D3 method of Grimme (available as of VASP.5.3.4) IVDW = 12 DFT-D3 method with Becke-Jonson damping (available as of VASP.5.3.4) IVDW = 2|20 Tkatchenko-Scheffler method (available as of VASP.5.3.3) IVDW = 21 Tkatchenko-Scheffler method with iterative Hirshfeld partitioning (available as of VASP.5.3.5) IVDW = 202 Many-body dispersion energy method (MBD@rSC) (available as of VASP.5.4.1) IVDW = 4 dDsC dispersion correction method (available as of VASP.5.4.1) All methods listed above add vdW correction to potential energy, interatomic forces, as well as stress tensor and hence simulations such as atomic and lattice relaxations, molecular dynamics, and vibrational analysis (via finite differences) can be performed. Note, however, that these correction schemes are currently not available for calculations based on density functional perturbation theory. N.B.: The parameter LVDW used in previous versions of VASP (5.2.11 and later) to activate DFT-D2 method is now obsolete. If LVDW=.TRUE. is defined, IVDW is automatically set to 1 (unless IVDW is specified in INCAR).
侯老师说[2]:
1,至于选择哪种,没办法推荐。大部分都是与所研究的体系有关。或者自己以所研究的体系中一个特定例子(即一个研究的问题)测试不同修正的方法,然后选择那个最好的,再计算其他的。
看来我们要把这些参数到底咋回事都要搞清楚。读读原始文献了。
[1] http://muchong.com/t-11094344-1
[2] http://muchong.com/t-11111598-1
继续搬砖。查找网络资源发现,Sobereva当年也仔细研究过这个。
DFT-D色散校正的使用[3]
文/Sobereva @北京科音
First release: 2013-Nov-7 Last update: 2018-Jun-5
多数传统的交换相关泛函,诸如B3LYP,由于相关势的长程行为不对,完全不能描述色散作用。因此用于研究色散作用主导的问题结果差得一塌糊涂,比如物理吸附、长链烷烃之类的大分子构象、弱极性分子团簇等等。解决密度泛函对色散作用描述能力很差的最有效的方法就是引入经验的色散校正项。不同的人提出过不同的色散校正方法,最成功也是目前最为流行的是Grimme提出的DFT-D。DFT-D实现非常容易,已在量化程序中得到广泛支持,它几乎不令计算量有任何增加。DFT-D有三个版本,目前最新、最严谨、参数最全(支持从H到Pu)的版本是在JCP,132,154104(2010)中提出的DFT-D3;老一些的DFT-D2是在JCC,27,1787(2006)中提出的,实际上性能与DFT-D3相仿佛,但参数少得多。更早的则是JCC,25,1463(2004)中提出的DFT-D1。
DFT-D对色散作用描述的改进立竿见影,例如笔者在J. Mol. Model., 19, 5387 (2013) http://link.springer.com/article/10.1007%2Fs00894-013-2034-2一文中比较了B3LYP加和不加DFT-D3校正时对氢气、氮气二聚体的计算结果,数据表明B3LYP根本没法得到二聚体稳定构型,或者说实际存在的构型下算得的二聚体相互作用能都是正的。而加上DFT-D3修正后,即B3LYP-D3,对这些二聚体的相互作用能计算结果则与金标准CCSD(T)起码定性一致。即便是静电主导的弱相互作用,诸如氢键、卤键,尽管传统泛函也能勉强凑合用,经过DFT-D校正后也能令计算精度显著提高。由于DFT-D同时也改进了中程的相关作用的描述,因此顺便对DFT泛函的热力学性质的计算精度也带来少许改进。目前来看可以说加入DFT-D校正有益无害,建议总加上DFT-D。自从DFT-D近几年流行起来后,曾经密度泛函对弱相互作用描述不好的老黄历彻底被颠覆,密度泛函摇身一变成了最为有效的计算弱相互作用的方法,特别针对是较大体系。同时大量相关的测试文章的涌现也让很多人清楚认识到了传统泛函对弱相互作用有多差。
原则上说,DFT-D3可以与任何交换相关泛函相结合。甚至对于弱相互作用描述已经很不错的双杂化泛函、M06-2X,加上DFT-D3后性能也能得到稍微的改进。但是也有一些出名的泛函本身就标配了特定形式的色散修正项,比如B97D、ωB97XD和B2PLYPD,在泛函定义的时候就专门给它们标配了DFT-D2形式的校正,显然就不能再给它们加上DFT-D3修正了。泛函结合DFT-D3校正后通常用“泛函名”+“-D3”来称呼,例如BLYP结合DFT-D3就叫BLYP-D3。
关于DFT-D在笔者另外的帖子里有更多讨论,见《乱谈DFT-D》(http://sobereva.com/83)和《谈谈“计算时是否需要加DFT-D3色散校正?”》(http://sobereva.com/413)。本文的重点是介绍一下DFT-D校正能如何通过Grimme的DFT-D3程序来计算,以及DFT-D如何在Gaussian、ORCA和GAMESS-US中启用。但在此之前先把DFT-D3的形式介绍一下。
DFT-D3实际上有两个版本,差异在于阻尼函数形式。阻尼函数用来调节色散校正在进程、中程距离时的行为,以避免double-counting问题(传统DFT泛函能够较好描述近距离的相关作用,如果近距离时校正能还较大的话就导致重复了)。DFT-D3原文当中用的是零阻尼(zero-damping)形式,这也是通常说的DFT-D3。而后来在JCC,32,1456(2011)中通过比较,发现使用物理意义更明确的Becke-Johnson阻尼(BJ-damping)可以让结果稍微更好点,对分子内色散作用的描述优势更显著些,这种校正形式通常被称为DFT-D3(BJ)。不过多数文章中在使用DFT-D3校正时并不做显著区分,表面写的是DFT-D3但实际上可能用的是DFT-D3(BJ)。
DFT-D校正能加到原先泛函计算的体系能量上就是校正后的能量。基于零阻尼的DFT-D3校正能写为
其中R_AB代表AB原子间的距离,上标n就代表了距离的n次方。C是原子间色散校正系数,依据一定规则进行计算。s_n是刻度因子。零阻尼函数f的表达式为:
其中R0_AB是原子对儿的截断半径,定义为√(C8_AB / C6_AB)。sr,n是刻度因子,γ是预设常数。之所以叫零阻尼,是因为随原子间距离R_AB减小阻尼函数逐渐衰减为0,使得DFT-D3校正能在较近距离时精确为0。
零阻尼形式对于每个泛函有4个可调参数,即s6、s8、sr,6、sr,8。实际上所有泛函的sr,8都是1。对于普通泛函,s6也为1。所以DFT-D3用于普通泛函时只需要拟合两个参数s8和sr,6即可。对于双杂化泛函,s6是一个小于1的值,也需要进行拟合来确定。
基于BJ阻尼的DFT-D校正能的形式为
对于每个泛函有四个参数,即a1、a2、s6、s8。依然是对于普通泛函s6为1,故一般只需要拟合三个参数,而双杂化泛函还得拟合它。
BJ阻尼使得中、近程距离的色散校正能虽很小,却并不接近于0。对于HF也可以用DFT-D3校正(HF对色散作用能完全无力),但只能用BJ而不能用零阻尼形式。对于明尼苏达系列泛函,如M06-2X,只能用零阻尼而不能用BJ阻尼,因为这类泛函已经表现出了中程相关,所以要用零阻尼形式来避免因引入中程校正而导致的double-counting。
上面的式子表明DFT-D3校正是基于原子对儿的,每一对儿原子对校正能有各自的贡献,但实际上的色散作用不是精确的对儿可加和的,还有多体项。于是DFT-D3还包含了三体项对校正能贡献。但由于对结果影响很小,所以考虑DFT-D3校正时一般不计算三体项。
由于DFT-D的形式简单,一阶和二阶解析导数很容易得到,所以原先DFT能做到几阶解析导数在校正后依然能做到几阶。不过如果把三体项考虑进去就得不到解析导数了。
Grimme的DFT-D3程序是用来计算DFT-D色散校正能的,程序一直在更新,下载地址见https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3。用户只需要提供原子坐标,程序就能立刻算出校正能以及校正能的导数。对于常用的泛函,DFT-D3的参数基本上都已经齐全。
编译过程很简单。对于Linux系统,下载后解压,进入目录输入make命令即得到了dftd3可执行程序,默认调用的是ifort来编译。对于Windows系统,假设使用的是Intel visual fortran,编译方式是新建项目,然后把压缩包内的.f文件以及param文件都加入项目然后编译即可。笔者编译好的Linux 64bit版和Windows版的DFT-D3程序都可以在这里下载:/usr/uploads/file/20150609/20150609203011_54761.rar。压缩包内带有.exe后缀的是windows版,不带后缀的是Linux版。
程序压缩包里包含手册,对运行参数介绍得很清楚。在Linux下,下面命令计算的是DFT-D3(BJ)对test.xyz里的体系的针对BLYP泛函的色散校正能
./dftd3 test.xyz -func b-lyp -bj
如果你用的是Windows,需要进入DOS然后同样以如上方式通过命令行指令来使用DFT-D3。因为DFT-D3程序没有交互式的界面,因此无法靠双击dftd3.exe图标来运行,必须通过命令行来运行!
程序会输出一堆参数,最后给出校正能
Edisp /kcal,au: -3.3486 -0.00533628
以a.u.为单位的色散校正能会同时写入到当前目录下的.EDISP文件中。
例子中的test.xyz是体系的xyz格式的坐标:
6 ltwd C 0.00000000 0.42021400 0.00000000 H -0.45599500 1.43114200 0.00000000 O 1.20173900 0.23372100 0.00000000 N -0.94137600 -0.56316900 0.00000000 H -0.64194700 -1.52713700 0.00000000 H -1.92633700 -0.35287300 0.00000000
第一行是原子数,第二行是标题,之后就是原子名和以埃为单位的坐标。
-func后面接的是所用的泛函,必须按照turbomole的泛函写法来写,诸如b3-lyp、pbe0、pbe、b2-plyp、cam-b3lyp、b2gp-plyp、m062x。如果不知道怎么写,可以到dftd3.f源文件里去搜b3-lyp,然后就能看到各种泛函的写法以及相应的参数,比如会看到
... case ("tpss0") rs6 =0.3768 s18 =1.2576 rs18=4.5865 case ("pbe0") rs6 =0.4145 s18 =1.2177 rs18=4.8593
...略
-bj代表用BJ阻尼,-zero代表用零阻尼。若写-old代表用DFT-D2形式。
DFT-D3程序默认不计算三体项,如果要计算的话就加上-abc,这对总校正能影响实在甚微,一般没必要计算。
加上-grad的话程序还会计算各个原子的色散校正能的梯度并写入到当前目录下的dftd3_gradient文件中。
加上-anal的话程序还会输出每个原子对儿对色散校正能的贡献。例如
analysis of pair-wise terms (in kcal/mol)
pair atoms C6 C8 E6 E8 Edisp
1 2 6 1 0.889542D+01 0.166327D+03 -0.1088 -0.1480 -0.25679
1 3 6 8 0.177185D+02 0.428060D+03 -0.1695 -0.2748 -0.44422
1 4 6 7 0.198682D+02 0.501857D+03 -0.1813 -0.3033 -0.48467
...略
就说明比如1号和3号原子之间的作用对色散校正能的贡献是-0.44422 kcal/mol。
用-anal的时候可以定义片段。方法是在当前目录下写一个名为fragment的文件,每行代表每个片段包含的原子,例如内容为
2-4
1,6
这代表第一个片段包含2、3、4号原子,第二个片段包含1、6号原子。剩下的第5号原子就自动作为第三个片段。程序会输出
group # atoms
1 2- 4
2 1 6
3 5
group i j Edisp
1 -- 1 -0.74
2 -- 1 -1.58
2 -- 2 -0.25
3 -- 1 -0.43
3 -- 2 -0.35
3 -- 3 0.00
由此可见片段1内的3个原子彼此间色散校正能是-0.74 kcal/mol,1、2号片段之间色散校正能是-1.58 kcal/mol。片段3就一个原子,所以片段内的没有色散校正能。以上数值加和就是总的校正能。
对于诸如B3LYP这样的泛函,由于完全表现不了色散作用,因此对于诸如B3LYP这样的泛函的DFT-D3校正能,特别是零阻尼下的,就可以近似视为是色散作用能。因此,利用DFT-D3程序我们可以直接计算出分子间色散作用能的大小。例如我们计算水二聚体,将两个水分别定义成两个片段,那么DFT-D3程序在-func b3-lyp -zero时给出的这两个片段间的校正能就可以当成是这两个水之间的色散作用能。尽管这种做法不怎么严格,Grimme也并不鼓励这么做,但从实际结果上看这么做基本还算合理。可以算是个计算分子间,或者两个分子当中的任意指定的片段之间色散作用能的便利的方法,但精度只能算得上是定性级别的。想准确计算分子间色散作用能,最可靠的方法是做SAPT(对称匹配微扰理论)得到相互作用能的色散成份,但这一方面比较耗时、麻烦,另一方面也没法简单分离出指定局部片段间的色散作用能
谈谈"计算时是否需要加DFT-D3色散校正?"[4]
笔者经常被问及的一个问题就是“我要计算xxx,是否需要加DFT-D色散校正?”。关于DFT-D,在以下博文里已经有充分的介绍:
乱谈DFT-D
http://sobereva.com/83
DFT-D色散校正的使用
http://sobereva.com/210
大体系弱相互作用计算的解决之道
http://sobereva.com/214
其实把以上博文所述内容都弄明白,什么时候该加DFT-D色散校正就已经应该心里很有数了。但是鉴于还是有大量初学者犯糊涂,我觉得有必要在此文里再强调一下,也算是对上面文章的补充。下文第一节是介绍一些常识,如果想对什么时候该加DFT-D这个问题了解更深刻可以看看,如果只是急于获得答案,可以略过。
化学体系中的相互作用主要分为化学键作用和弱相互作用两大类。从相互作用能上看,前者一般强度比较强,共价键和离子键都属于此类,也可以叫强相互作用;而后者一般强度比较弱,比前者通常弱一个数量级。弱相互作用既可以是分子间的,也可以是分子内的。弱相互包括范德华作用、pi-pi堆积作用、氢键、二氢键、卤键,以及后来很多人炒作概念而提出来的碳/硫/磷/金/银/铜...键等等。还有个词叫做非共价相互作用(noncovalent interaction),这个词的范畴相当于弱相互作用和离子键的并集。
弱相互作用如上所列,形式多样,但主要本质只有两个:
(1)静电相互作用。可以起到互斥作用也可以起到吸引作用,看具体情况(体系,以及相对位置)。这里说的静电相互作用也把极化作用包含进去了。
(2)色散作用。它起到吸引作用。必须从量子力学角度才能予以解释。而从量化理论角度来讲,它对应于电子的长程的库仑相关作用。
一般强度的氢/二氢/卤/硫/磷键等等是以静电吸引主导,色散作用为辅构成的,这点通过能量分解也可以体现出来。而一般说的范德华相互作用,以及pi-pi堆积作用的本质都是色散作用。
注1:还有个作用叫交换互斥,它起到互斥作用,不管算什么类型相互作用都要考虑这个,只有在较近距离(比如小于两个原子的范德华半径和)的时候才开始凸显出来,且原子间距离越近此作用越大。正是这个作用使得弱相互作用的势能曲线总是有极小点,而不会令原子间距离因为静电吸引和色散作用而无限减小。
注2:“范德华作用”这个词描述的是电中性的两个原子之间的非共价相互作用,故本质包含了交换互斥和色散吸引作用两个部分。这个词可以囊括比如两个氩原子间的相互作用的全部内涵。这个词也可以用来描述无极性分子间的全部相互作用(但实际上,没有绝对严格意义的无极性分子,比如就连一般说的无极性分子氮气也有四极矩,可以认为是局部极性,因此氮气之间也有静电相互作用,讨论见http://sobereva.com/209。但本文就不咬文嚼字了)。至于一些初级的教材里说的范德华作用包括“诱导力、色散力、取向力”,这种歪曲的说法应该彻底从量化研究者的脑海中删除。
注3:弱相互作用实际上也伴随着一些共价作用,但共价作用此时仅仅是个“陪衬”而已,对弱相互作用能贡献程度极小,故在弱相互作用范畴里不需要专门特意考虑(但是弱相互作用范畴里面也有强度很大的情况,比如能到100kJ/mol数量级的氢键,此时共价作用就不能忽视了,此文说的弱相互作用不考虑此情况)。
不管是什么DFT泛函,什么理论方法,对上面提到的静电相互作用和交换互斥作用总是可以至少定性正确地表现。然而老一代泛函对于色散作用描述极糟(如SVWN、PBE、PW91),或者根本不能描述(如BLYP、B3LYP),因此这些泛函完全没法描述范德华结合、pi-pi堆积这样本质来自于色散吸引的弱相互作用,结果会定性错误;而考察氢键等静电主导的弱相互作用时,虽然多数情况定性正确,但肯定定量明显不准。这些泛函有这些问题一方面是拟合参数的时候没有考虑弱相互作用体系,一方面是泛函的相关势的长程渐进行为明显不对。
为什么Grimme要提出DFT-D?出发点就是解决老一代泛函计算弱相互作用糟糕的情况。更确切来说,解决的是对色散作用描述糟糕的问题。因此,DFT-D叫做色散校正。但要注意,色散校正,或者解决传统泛函描述色散作用烂的方法绝对不止DFT-D这一种,还有很多其它做法,比如vdW-DF、XDM、TS、LRD、DFT-ulg、dDSC、OBS等等,由于这些方法要么实现复杂、要么会导致计算耗时增加许多、要么效果不好、要么还没有被主流程序普遍支持等原因,时下流行程度和DFT-D相比微不足道,以至于DFT-D简直成了色散校正的代名词。
DFT泛函不加色散校正就不能描述好弱相互作用了?非也。由于弱相互作用越来越受重视,从21世纪00年代中期开始提出的泛函越来越多地注重对弱相互作用的描述。比如2006年的M05-2X描述各种弱相互作用就已经挺不错了,因为在拟合参数的训练集中已经纳入了弱相互作用体系,之后08年的M06-2X在这方面又得到了一定进步。而比如08年提出的ωB97XD,在提出当初就已经带了DFT-D2校正以使得计算各类弱相互作用能够有比较好的能力。至于从2006年开始提出的双杂化泛函,由于里面已经包含了MP2成份,而MP2本身就已经能够定性正确描述色散作用,因此双杂化泛函对各类弱相互作用都有比较好的描述能力,哪怕其中有的不及M06-2X,但也总比较老一代的泛函要强得多得多。
DFT-D第一代是2004年提出的,最新的是第四代(DFT-D4)。时下,最火的还是2010年提出的DFT-D3,它纯粹基于几何结构来计算色散校正能,加到原有泛函的能量上就是色散校正后的能量,比如加到B3LYP能量上就成了B3LYP-D3能量。更多信息见前面提到的《DFT-D色散校正的使用》和《乱谈DFT-D》,这里不再累述。DFT-D1/2/3都只利用几何信息而不考虑体系电子结构,计算都很迅速,几乎不增加任何计算量,故是“免费”的。而到了DFT-D4,校正能才考虑电子结构(依赖于原子电荷),计算起来更为复杂,也不再是完全“免费”的了,但带来的好处就是可以充分体现电子结构对色散校正的影响(比如同一个结构下,激发态和基态是不同的,单重态和三重态是不同的,中性和离子状态是不同的...特别是对于计算诸如金属氧化物等原子电荷数值较大的情况,由于考虑了实际电子结构,结果会比用DFT-D3更好)。
对于该不该、有无必要加DFT-D3这个问题,下面对计算弱相互作用的情况和计算其它任务的情况分别说。
这里说的计算弱相互作用,显然包括计算分子间相互作用能、优化分子二聚体/多聚体几何结构(以及随后的振动分析)、研究物理吸附等场合。对于很多柔性大分子,存在显著的分子内相互作用,因此弱相互作用对其不同构象间的相对能量、构象结构影响极大,故研究它们的构象问题也属于下述的“计算弱相互作用”的范畴里。是否存在明显分子内/分子间相互作用其实一看结构心里就有数,只要非键连的原子之间距离挨得较近,就应当考虑弱相互作用。如果化学直觉不够而不好判断,可以用Multiwfn做IGM分析,见《通过独立梯度模型(IGM)考察分子间弱相互作用》(http://sobereva.com/407);或者做RDG/NCI分析,见《使用Multiwfn图形化研究弱相互作用》(http://sobereva.com/68),操作都十分简单,敲几下键盘的事。
当你研究的问题涉及到弱相互作用,而且当前用的理论方法不足以合理描述色散作用,色散校正是一定要加的!鉴于DFT-D3是免费的而且是最流行的色散校正方法,我们一般首选加DFT-D3。比如说B3LYP描述色散作用完全失败,当拿它计算pi-pi堆积这种完全由色散作用所控制的弱相互作用时,不加色散校正结果是荒唐的。而用B3LYP计算氢键这种静电主导而色散作用为辅的弱相互作用时,加色散校正后的结果也才可能用来发文章(可能有人看到这句话会突然跳出来说:我看那谁谁谁,还挺有名的,以前用B3LYP算氢键不是也发文章了?要知道,弱相互作用的计算在最近十几年发展极其迅速,05年之前DFT算弱相互作用的文章完全没法看,05~10年的文章也只能有选择性地参考,10年之后的DFT算弱相互作用的文章才算是进入成熟期。研究者们对DFT做弱相互作用计算的认识在这期间也发生了翻天覆地的变化。除非碰上外行审稿人或者投极烂的期刊,这年头用B3LYP算弱相互作用能发表才怪。)
当研究弱相互作用相关的问题时,如果目前用的DFT泛函已经可以不错地描述色散作用了,比如用的是M06-2X,那显然DFT-D3不是必须加的,加了不会有什么坏处,但带来的改进甚微。如果用的是双杂化泛函,由于对色散作用也已经描述不错了,因此不加D3算各种弱相互作用的精度也不赖,至少不至于因不加D3而被审稿人拒,但考虑到从统计结果来看,加了D3后精度改进还是有不少油水的,所以建议总是加上。很多后来提出的双杂化泛函,比如PWPB95、DSD-PBEP86,由于在提出时就已经带着D3了,所以用的时候更是强烈建议带着D3(值得一提的是,Gaussian16里写DSD-PBEP86关键词的时候自动就是带DFT-D3的)。
以下给出Grimme在Phys. Chem. Chem. Phys., 19, 32184 (2017)中列出的几个泛函的弱相互作用平均计算误差,D3用的是BJ阻尼:
BLYP:18.56 kcal/mol
BLYP-D3:3.70 kcal/mol
B3LYP:15.43 kcal/mol
B3LYP-D3:2.94 kcal/mol
PWPB95:5.74 kcal/mol
PWPB95-D3:2.27 kcal/mol
可见,BLYP、B3LYP加D3后改进甚大,不加D3误差大到根本没法用。PWPB95双杂化泛函加了D3之后也有很大改进。
下面是不同明尼苏达系列泛函加D3前后的平均误差情况。明尼苏达系列泛函只能用零阻尼形式的D3校正。
可见,对于计算弱相互作用常用的M06-2X,加D3的意义不大,改进甚微。仅对于M06-2X计算pi-pi堆积,我建议还是尽量加上D3,因为对这种情况M06-2X有时误差可能略大,加D3后有不可忽略的改进。并非所有明尼苏达系列泛函算弱相互作用都不错,有的实在是非常糟糕,而且加D3都救不活,比如MN15L这种。关于某些明尼苏达系列泛函加D3校正后结果仍较烂原因的一些细致的探讨,可以看J. Phys. Chem. Lett., 6, 3891 (2015)。
对于后HF、多参考方法等等,是没有色散校正这一说的。DFT-D仅能用于DFT、DFTB、HF、半经验。高档次后HF方法,比如MP4、CCSD(T),已经可以很准确计算各类弱相互作用体系了,也没什么系统性误差,且误差已明显小于DFT泛函结合DFT-D校正后的误差了,显然对这些精细的方法再试图通过相对来说形式较为粗俗的DFT-D来校正是没有任何意义的,因此也没人对这些方法去拟合DFT-D参数。而至于MP2,由于耗时和双杂化差不多,而算弱相互作用精度又不及双杂化(氢键除外),因此也没人专给它搞DFT-D参数。
有人问:“我用xxx泛函,我知道计算弱相互作用能的时候有必要加D3,那优化和振动分析有没有必要?”显然有必要!D3校正能是依赖于几何结构的,因此,D3校正本质上校正的是体系的势能面,会令所有依赖于势能面的问题都受到影响,自然也包括优化出的结构(势能面极小点的位置)、振动频率(取决于势能面极小点处Hessian矩阵)、IRC(质权坐标下的势能面上的能量极小路径)等。你算能量时候用D3,就代表你已经知道考虑D3校正对于研究当前体系是有意义的,如果几何优化等任务不加D3,怎么说得通?必然会糟至审稿人吐槽,何况D3还是免费的。所以,优化、振动分析、单点、IRC等等,要加D3都加,如果确信没必要加,那就都不加。
另外,如果你研究的是一系列体系,用的是同一个泛函,只要其中一个加了D3,其它所有的计算全都得加D3,否则没法横向比较。如果涉及到能量直接求差的情况,比如复合物能量减去各个片段能量、一个构象减去另一个构象的能量,那所有涉及的体系加不加D3就更得统一,否则能量求差毫无意义。
可能有时候研究的问题并不直接牵扯弱相互作用,但是体系内又存在弱相互作用,这时候如果你用的泛函描述色散能力不足,那照样也得加D3。比如说研究某个大分子的局部异构化过程,这个过程本身不怎么涉及到弱相互作用特征的变化,但大分子其余部分存在分子内弱相互作用(比如有一些内氢键),对这种情况,你若是用M06-2X算的则D3可加可不加,但用B3LYP算就应当加D3。或许你不觉得对此问题加D3有什么必要,但毕竟加D3又不花钱又基本没坏处,还省得之后审稿人找碴又让你加D3重算,何故不加?而且还完全排除了不考虑色散作用会对结果产生不可忽略的误差的风险。
因为D3校正仅仅是基于几何结构计算校正能,与电子态无关,而且只改变能量而不影响体系的波函数,因此凡是只涉及单个几何结构的任务都不需要加D3,或者说,加不加结果都一样。比如计算某个结构下的gap、轨道、偶极矩、极化率、NMR、原子电荷、键级、静电势分布,以及计算垂直激发能、理论电子光谱(不考虑振动耦合时)、垂直电离能、垂直电子亲合能、垂直单-三重态能量差等,D3丝毫不影响结果。但是,D3可以间接地对它们产生影响。比如一个二聚体,用B3LYP和B3LYP-D3优化出来的结构可能相差很大,那显然最终计算的偶极矩、吸收光谱等性质也会有明显不同。
诸如有人问过我,计算激发态用不用加D3,那显然得看你算激发态的什么。算垂直吸收/发射能和振子强度、计算激发态与基态的密度差、计算NTO等等,前面说了,D3丝毫不可能影响结果,因为结构没变,而且D3不影响波函数。而对于激发态几何优化,D3显然会影响结果,而且如果优化基态时加了D3则优化激发态也得加,得保持统一。
对于主要只体现化学键作用的体系,加D3是没必要的,虽然加了也没坏处。诸如优化有机小分子、简单配合物(如二茂铁)、简单阴阳离子体系(比如乙酸钠、硫酸钙)等等,加D3基本不会对结果产生任何影响。一方面是D3里面有阻尼函数,阻尼函数会使得校正函数在近程区域衰减为0,从而避免影响原先DFT泛函就已经能描述良好的范畴,另一方面是对于这些较强的相互作用,微小的色散作用的影响本身就是可以忽略的。但如果是有机大分子、配体比较大的配合物(此时配体间容易有弱相互作用)、比较大的阴阳离子对等等情况,弱相互作用的存在不可忽略,当使用对色散作用描述烂的泛函的时候建议加上D3。
对于计算小体系化学反应能(也包括原子化能、生成焓等),没有太大必要加D3,但是如果体系比较大,对于描述色散作用能力不足的泛函还是建议加上D3。毕竟体系原子数一多,体系内弱相互作用总量就会较大,更严谨地来说,体系中的中、长程电子相关总效应会较大。D3虽然原本目的是校正色散作用的描述能力,但实际上也等效地一定程度改进了很多老泛函在中、长程区域对电子相关描述能力的不足。
下面的图来自Phys. Chem. Chem. Phys., 19, 32184 (2017),可见D3对于反应能的计算精度是有一定改进的,而且对大体系改进得比小体系更多。
但加D3后结果变差也不是不可能,要想找反例总能找得到。比如JCTC, 13, 3537 (2017)报道PBE0-D3算生成焓比PBE0精度有所下降。不过,相对于一些负面消息,更多的文献测试还是体现出D3会带来积极的效果。
从上图看D3对能垒计算结果并无改进,但实际研究工作往往是研究一整套化学反应过程,若优化极小点、算极小点能量用D3的话,对过渡态也一定得用D3来保持统一。
虽然上面已经讲了很多,但肯定有些初学者由于基础知识不足,面对一些情况对于加不加D3还是会迷糊。简单一句话:不知道该不该加就加。反正又免费用着又方便,绝大多数情况都有益无害。
持续更新中。。。。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-12-28 13:10
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社