用NWChem做SODFT在DFT计算中考虑旋轨耦合效应
用NWChem做SODFT在DFT计算中考虑旋轨耦合效应
Use NWChem to do SODFT and consider spin-orbit coupling effect in DFT calculations
文/Sobereva @北京科音 2017-Apr-12
Use NWChem to do SODFT and consider spin-orbit coupling effect in DFT calculations
文/Sobereva @北京科音 2017-Apr-12
1 前言
众所周知,对重原子构成的体系,尤其是从第五周期开始(H,He算第一周期,下同),需要考虑相对论效应,否则结果可能严重错误。相对论效应分为标量和旋轨耦合两部分。标量相对论效应好考虑,用相对论赝势时即可等效体现,多数主流量化程序如Gaussian也支持标量相对论哈密顿,结合专为标量相对论计算优化的全电子基组即可(见《在赝势下做波函数分析的一些说明》http://sobereva.com/156)。比较棘手的是旋轨耦合的考虑,本文只说单个态的旋轨耦合问题。最直接的考虑方式是做二分量相对论计算,但比较耗时,且支持它的ADF、Turbomole都是收费的,而免费的Dirac等程序则比较小众且功能局限性大。Gaussian里有个DKHSO方法可以在单点计算时考虑旋轨耦合效应,但是只能用于HF,因此没什么实用性。一种常用的又便宜又好的考虑旋轨耦合的做法是使用相对论赝势,但不仅要包含我们平时常用的标量势,还同时要考虑旋轨势,这种方式做DFT称为Spin-orbit DFT (SODFT),Molpro、Turbomole、NWChem都支持。本文介绍使用其中免费的NWChem程序做SODFT的方法,并会以TlCl和HgI的解离能计算做为例子。在NWChem中,SODFT有解析梯度,因此能够容易地优化分子结构,不过Hessian只有数值的(只考虑标量势的话倒是有解析Hessian)。NWChem的SODFT还一个不足之处是没法利用对称性加速计算。
NWChem的安装方法看此文:《NWChem 6.6编译方法》(http://sobereva.com/270)。本文使用NWChem 6.6版。
2 SODFT计算时赝势和赝势基组的选择
赝势包括标量势、旋轨势和核极化势。标量势部分是我们最常使用的,所有相对论赝势都提供了标量势,但不是所有赝势都给出了旋轨势。给出了旋轨势的赝势比较常见的有Stuttgart赝势和CRENB赝势,对周期表覆盖得都比较全面。后者很老,精度也不比前者好,所以本文将会用前者。Stuttgart赝势也分多种,包括HF系列(赝势不体现相对论效应),WB系列(基于Wood-Boring准相对论计算的数据搞的)和DF系列(基于Dirac-Fock全相对论计算的数据搞的),我们一般都用M版(拟合赝势时考虑多个价电子),如MDF。除赝势外还要考虑赝势基组。基于Stuttgart赝势的常用的几种选择有:
1 Stuttgart赝势的标配基组。缺点是对于大多数元素,不像下面几种给出了不同档次基组的定义,选择余地小。
2 cc-pVnZ-PP系列:结合的是Stuttgart小核MDF赝势
3 def2系列:从第五周期开始是赝势基组,结合的是Stuttgart小核赝势,对s、d、镧系元素用的是MWB,对p元素用的是MDF
4 dhf系列:和def2都是Weigend主要搞的,只定义了五、六周期的s和d族元素,p族元素和def2定义相同。此系列分为dhf和dhf-2c两种,前者主要用于做一分量计算(即常规的赝势下的计算),把def2对s、d元素用的标量势从MWB改为MDF并重新调整了赝势基组定义,由此使结果略好了一点点;后者在前者基础上增加了额外的p,d基函数,使得带着旋轨势做计算时结果更好
以上赝势基组中,做SODFT时能用dhf-2c就用dhf-2c,毕竟是给此目的专门优化的。对于本文要考察的元素,Tl、Hg、I都可以用dhf-2c,我们用其中TZVP级别的就够了,即dhf-TZVP-2c。虽然也有更好的dhf-QZVP-2c,但对于DFT计算升到QZ级别意义并不大。对于Cl,由于相对比较轻,其引起的旋轨耦合效应不会给结果带来明显影响,所以我们不对它用赝势,而对它用def2-TZVP全电子基组。
3 计算TlCl
本文我们要尝试重复下面图中的TlCl和HgI的数据,数据来自Theor. Chem. Acc., 130, 633-644 (2011)。图中BDE是键解离能,单位是kcal/mol。键解离能的计算是用解离后的原子的焓减去分子的焓。图中BDE下面的列是文中只考虑标量相对论的计算结果,BDE+SO是考虑旋轨耦合后的。![](/usr/uploads/image/20170412/20170412235424_10703.png)
我们先算TlCl。由于文中用的是PBE0,所以我们也用PBE0,此泛函对这类体系是比较合适的。我们用的键长直接取实验值,大量双原子分子的实验平衡距离可以在NIST网站http://webbook.nist.gov/chemistry/上很容易地查到,注意勾选Constants of diatomic molecules,之后数据页面里r_e就是平衡距离了,对此分子是2.484826埃。
首先我们得获取赝势和赝势基组定义。def2和dhf系列在Turbomole基组库http://www.cosmologic-services.de/basis-sets/basissets.php里是最全的(EMSL上的不仅不全,而且没有旋轨势)。在此网站里面选上Tl,然后选dhf-TZVP-2c,输出格式就用默认的Turbomole(没直接提供NWChem的比较可惜)。给出的数据包括基组部分、标量势部分和旋轨势部分。
PS1:值得一提的是,当选择dhf-TZVP-2c的时候,对于非第五、六周期的元素,自动给出的是def2对应档次基组的定义
PS2:此网站输出格式里只有选Turbomole或Molpro这两种支持旋轨势的程序时输出内容才有旋轨势,我们这里不选Molpro是因为没有Turbomole的格式那么容易手动改成NWChem的格式
然后我们把Turbomole格式的基组、标量势和旋轨势的定义都改成NWChem的格式,不知道NWChem的格式怎么写就RTFM。一个较大区别是Turbomole里把收缩系数放在第一列,指数放在最后一列,和NWChem的顺序是反过来的。最终TlCl的SODFT任务的输入文件如下
start
PRINT low
geometry noautosym
Tl 0.00000 0.00000 0.00000
Cl 0.00000 0.00000 2.484826
end
basis "ao basis" spherical
Cl library def2-TZVP
Tl s
729.65038145 0.13672829828E-03
46.665548707 0.60443439951E-02
20.970448726 -0.20022066697
14.149588677 0.40801678488
Tl s
20.730134285 -0.71861135918E-01
6.1527631309 0.98057508445
Tl s
1.5757324480 1.0000000000
Tl s
0.74980169458 1.0000000000
Tl s
0.19536816194 1.0000000000
Tl s
0.70878767298E-01 1.0000000000
Tl p
15.383852616 0.61717949180
14.814929544 -0.72859235151
6.7261253658 0.40438195364
Tl p
1.9626182149 0.43157661160
1.0331857851 0.39230403853
0.53837445996 0.14007406820
Tl p
0.24446611676 1.0000000000
Tl p
0.90785377260E-01 1.0000000000
Tl p
0.33401321498E-01 1.0000000000
Tl d
57.606819928 0.16054811114E-03
9.7368866667 0.24456562496E-01
6.9256201679 -0.69914775031E-01
2.1396230731 0.19496269490
1.0836187110 0.29731629705
0.52356298209 0.23728708102
Tl d
0.23842309375 1.0000000000
Tl d
0.95000000000E-01 1.0000000000
Tl f
0.28435 1.0000000000
Tl f
0.95 1.0000000000
Tl p
32.0 -0.0005
16.0 -1.0
Tl p
8.0 -1.0125
4.0 0.7450
Tl p
2.0 1.0
Tl d
3.0 1.0
END
ECP
Tl nelec 60
Tl ul
2 1.00000000 0.00000000
Tl s
2 12.16780500 281.28466300
2 8.29490900 62.43425100
Tl p
2 7.15149200 4.63340800
2 5.17286500 9.34175600
2 9.89107200 72.29925300
2 9.00339100 144.55803700
Tl d
2 7.13021800 35.94303900
2 6.92690600 53.90959300
2 5.41757000 10.38193900
2 5.13868100 15.58382200
Tl f
2 5.62639900 15.82548800
2 5.54895200 21.10402100
2 2.87494600 2.91512700
2 2.82145100 3.89690300
Tl g
2 6.67905700 -7.49453400
2 6.70683500 -9.54057500
2 7.20928400 -7.79799200
2 7.07096400 -9.25952400
END
SO
Tl p
2 7.151492 -9.266817
2 5.172865 9.341756
2 9.891072 -144.598506
2 9.003391 144.558037
Tl d
2 7.130218 -35.943039
2 6.926906 35.939729
2 5.417570 -10.381939
2 5.138681 10.389215
Tl f
2 5.626399 -10.550326
2 5.548952 10.552010
2 2.874946 -1.943418
2 2.821451 1.948451
Tl g
2 6.679057 3.747267
2 6.706835 -3.816230
2 7.209284 3.898996
2 7.070964 -3.703809
END
dft
xc pbe0
end
task sodft
输入文件还是很容易理解的,照常定义基组部分,ECP段落指定标量势部分,SO字段指定旋轨势部分,最后task里写上SODFT就行了。noautosym写不写无所谓,由于SODFT不支持对称性,不写它时程序也不会利用对称性来加速。
然后我们运行此任务,输出文件里这部分即是SODFT能量:
Total DFT energy = -632.863447528456
在笔者的Intel 4核i7-2630QM机子下,做这个计算耗时是20.8s,若只考虑标量势是10.8s,可见考虑旋轨势还是会增加不少耗时,不过一般完全承受得起,总比做二分量相对论计算来考虑旋轨耦合效应便宜得多。
之后我们把TlCl的文件分别改写为Tl和Cl的文件,二者都是二重态,所以在dft字段里面要写上mult 2。算出的结果是Tl=-172.724034492809,Cl=-459.997396831618。因此基于单点能算的键解离能是627.51*(-459.997396831618-172.724034492809+632.863447528456)=89.116 kcal/mol。
到此还没完,我们还得计算焓的热校正量对解离能的修正。我们就不用SODFT算频率了,就直接在Gaussian里用# PBE1PBE/def2TZVP opt freq scale=0.982关键词对TlCl优化和做振动分析就完了,这里0.982是ZPE校正因子(为了省事我们不单独考虑升温对焓影响那部分的校正因子)。实际上我们查不到PBE0/def2-TZVP级别的校正因子,但由于J. Phys. Chem. A, 111, 11683 (2007)中给出PBE0结合6-311+G(2df,p)这样不错级别基组的校正因子是0.9824,所以有理由认为0.982对于PBE0结合def2-TZVP是基本合适的。计算后得到的焓的热校正量输出为
Thermal correction to Enthalpy= 0.004388
对于Tl和Cl没必要做振动分析,因为单个原子只有平动对焓有贡献,且是精确已知的,即5/2*RT,常温下对应0.002360 Hartree。所以常温下焓的热校正量对TlCl键解离能的贡献可以计算为627.51*(0.002360+0.002360-0.004388)=0.208 kcal/mol。
最后,得到TlCl的常温下键解离能是89.116+0.208=89.324 kcal/mol,这和前面贴出来的实验值88.1 kcal/mol相符极好!(对普通DFT泛函来说能达到1kcal/mol程度误差着实难能可贵)
4 计算HgI
HgI是二重态,解离后会生成单重态的Hg和二重态的I。HgI的键长在NIST上查不到,虽然某些文献里估计能查到,但这里不深究,就直接用PBE0/def2-TZVP级别来优化键长,之后再做单点计算。优化后键长是2.75663埃。(当然,若要求更精确,可以直接用SODFT优化。笔者也试了下,i7-2630QM上耗时336s,结果是2.76379埃。可见对此体系是否考虑旋轨耦合对结构影响很小,对TlCl测试过也是如此)。我们使用上一节的方式,获取Hg和I的dhf-TZVP-2c赝势基组定义以及对应的标量势和旋轨势,写成HgI的PBE0的SODFT任务的单点输入文件,如下所示:
start
PRINT low
geometry noautosym
Hg 0.00000 0.00000 0.00000
I 0.00000 0.00000 2.75663
end
basis "ao basis" spherical
I s
5899.5791533 0.24188269271E-03
898.54238765 0.15474041742E-02
200.37237912 0.42836684457E-02
31.418053840 -0.39417936275E-01
15.645987838 0.96086691992
I s
11.815741857 0.57815778954
6.4614458287 0.37374293124
I s
2.3838067579 1.0000000000
I s
1.1712089662 1.0000000000
I s
0.32115875757 1.0000000000
I s
0.12387919364 1.0000000000
I p
185.43362455 0.83127824000E-03
20.091408146 0.63991653000E-01
9.7577022390 -0.27791138000
I p
13.068307912 -0.49793790382E-01
3.5818714205 0.38212490511
2.0282441852 0.70447564804
1.0181492146 0.33781067803
I p
0.46673773115 1.0000000000
I p
0.19242597960 1.0000000000
I p
0.74508878495E-01 1.0000000000
I d
124.20341062 0.65671747209E-03
34.587311801 0.51648185674E-02
12.767328064 -0.19881371307E-01
4.7745100133 0.21386794109
2.4582209028 0.43405444707
1.1923708147 0.37850637882
I d
0.52883971906 1.0000000000
I d
0.17008164307 1.0000000000
I f
0.44141808 1.00000000
I f
2.18 1.0000000000
I p
12.0 -0.3
6.0 1.0
I d
2.0 1.0
Hg s
48.013786990 0.58613168385E-02
21.239875095 -0.17590367988
15.876100879 0.35780354753
Hg s
5.4837607070 1.0000000000
Hg s
1.5480592128 1.0000000000
Hg s
0.72425230437 1.0000000000
Hg s
0.16369906863 1.0000000000
Hg s
0.57211615334E-01 1.0000000000
Hg p
23.291760168 -0.83564430982E-02
13.028969569 0.83703058587E-01
6.5100040792 -0.31023833705
Hg p
1.8167935815 1.0000000000
Hg p
0.90079391013 1.0000000000
Hg p
0.41304090835 1.0000000000
Hg p
0.11845879331 1.0000000000
Hg p
0.36084184656E-01 1.0000000000
Hg d
15.176302343 0.60654575178E-02
6.7004896493 -0.59880306300E-01
1.9144256118 0.31411145584
0.88641552102 0.45081091161
Hg d
0.38364767725 1.0000000000
Hg d
0.14936891490 1.0000000000
Hg f
0.79569 1.0
Hg p
19.3582 0.0635
8.7992 -0.4607
Hg p
3.9996 0.8853
Hg d
2.7000000000 1.0000000000
END
ECP
I nelec 28
I ul
2 1.000000 0.000000
I s
2 40.033376 49.989649
2 17.300576 281.006556
2 8.851720 61.416739
I p
2 15.720141 67.416239
2 15.208222 134.807696
2 8.294186 14.566548
2 7.753949 28.968422
I d
2 13.817751 35.538756
2 13.587805 53.339759
2 6.947630 9.716466
2 6.960099 14.977500
I f
2 18.522950 -20.176618
2 18.251035 -26.088077
2 7.557901 -0.220434
2 7.597404 -0.221646
Hg nelec 60
Hg ul
2 1.0000000 0.0000000
Hg s
2 12.4130710 275.7747970
2 6.8979130 49.2678980
Hg p
2 11.3103200 80.5069840
2 10.2107730 161.0348240
2 5.9398040 9.0834160
2 5.0197550 18.3677730
Hg d
2 8.4078950 51.1372560
2 8.2140860 76.7074590
2 4.0126120 6.5618210
2 3.7953980 9.8180700
Hg f
2 3.2731060 9.4290010
2 3.2083210 12.4948560
Hg g
2 4.4852960 -6.3384140
2 4.5132000 -8.0998630
END
SO
I p
2 15.720141 -134.832478
2 15.208222 134.807696
2 8.294186 -29.133096
2 7.753949 28.968422
I d
2 13.817751 -35.538756
2 13.587805 35.559839
2 6.947630 -9.716466
2 6.960099 9.985000
I f
2 18.522950 13.451079
2 18.251035 -13.044039
2 7.557901 0.146956
2 7.597404 -0.110823
Hg p
2 11.31032000 -161.01396700
2 10.21077300 161.03482400
2 5.93980400 -18.16683200
2 5.01975500 18.36777300
Hg d
2 8.40789500 -51.13725600
2 8.21408600 51.13830600
2 4.01261200 -6.56182100
2 3.79539800 6.54538000
Hg f
2 3.27310600 -6.28600100
2 3.20832100 6.24742800
Hg g
2 4.48529600 3.16920700
2 4.51320000 -3.23994500
END
dft
mult 2
xc pbe0
end
task sodft
之后稍作修改,得到Hg和I的PBE0-SODFT任务输入文件。计算后将能量求差,得到基于电子能量的BDE:
627.51*(-295.746562477964-153.684090186717+449.448365018831)=11.115 kcal/mol
如前面的图所示,HgI的解离能的三个不同来源的实验值为7.8、8.1、8.9 kcal/mol,当前算的和哪个都想差甚远,因此PBE0对此体系不够给力。对这种情况,我们可以在DFT级别计算出旋轨耦合对BDE的校正,然后加到CCSD(T)/def2-QZVP这样只考虑标量相对论的高计算级别得到的解离能上。为此,我们把上面用的HgI、Hg、I输入文件里的sodft改成dft(之后SO段落可以删掉,不删也无所谓),这样计算就只利用标量势部分了,旋轨耦合效应就不考虑了。结果为
627.51*(-295.697077663054-153.500765444841+449.222674994391)=15.582 kcal/mol
因此,旋轨耦合对BDE的校正量为11.115-15.582=-4.467 kcal/mol。
在CCSD(T)/def2QZVP级别算的HgI解离过程的电子能量变化:
627.51*(-153.2231386-297.2662532+450.5085777)=12.039 kcal/mol
用Gaussian在PBE0/def2-TZVP下做振动分析,得到常温下焓的热校正量对HgI键解离能的贡献:
627.51*(0.002360+0.002360-0.004277)=0.278 kcal/mol
put together,得到HgI的常温下的BDE:
12.039-4.467+0.278 = 7.85 kcal/mol
可见结果和实验符合得极好!
5 SODFT方式计算BDE总结
前面HgI的例子很有代表性,说明想准确计算重原子的键解离能,最简单省事的过程就是:1 用较合适的泛函,在def2-TZVP或相似级别的结合相对论赝势的赝势基组(如cc-pVTZ-PP,但不划算)下进行优化和做振动分析,得到焓的热校正量对BDE的影响。注意应考虑ZPE校正因子
2 用较合适的泛函,在1的结构下用dhf-TZVP-2c赝势基组与配套的MDF标量势+旋轨势做SODFT计算,再去掉旋轨势做普通DFT计算,利用差值得到旋轨耦合对BDE的影响
3 基于1的结构,在尽可能高的级别下计算考虑了标量相对论效应的单点,由此得到基于电子能量算的BDE。体系很小时计算级别一般用CCSD(T)结合def2-QZVP或cc-pVQZ-PP赝势基组,赝势只考虑标量势部分。可以把结果再外推到CBS。
最后把1、2、3相加,就得到了特定温度下的键解离能。
上面说的是体系只含第五、六周期s,d,p族元素的情况。对于体系存在其它元素时,笔者认为应当这样考虑:
对前三周期的元素,就用全电子基组就行,不需要考虑相对论效应。
对第四周期元素,考虑相对论效应有益。分为两种情况:
(a)Cu之前的:Stuttgart赝势没有给出这些元素的旋轨势,cc-pVnZ-PP也没有给出这些元素的定义,def2对这些还是全电子基组。要考虑相对论效应就用Stuttgart赝势结合其标配赝势基组就完了。非要也考虑旋轨耦合那就用CRENB标量势+旋轨势结合CRENBL赝势基组吧。
(b)Cu~Kr:用cc-pVnZ-PP赝势基组,结合对应的MDF标量势和旋轨势。从EMSL上取cc-pVnZ-PP定义的时候赝势部分只有标量势部分,旋轨势得从Stuttgart赝势的官网http://www.tc.uni-koeln.de/PP/clickpse.en.html上自行取,输出格式必须选molpro,要选小核的MDF版本(比如Br要选ECP10MDF)。输出的赝势里前几行是标量势(和EMSL上cc-pVnZ-PP直接带的一致),后几行是旋轨势。
对于镧系、锕系及一些第七周期的元素,就用Stuttgart赝势官网上小核标量势+旋轨势结合标配的赝势基组。情况较复杂。有的有MWB有的有MDF有的都有(优先用后者),有的有旋轨势有的没有旋轨势(没有的话就没辙了)。
另外,上面说的第1步的优化+振动分析,如果求准确,可直接用SODFT来做,但是由于SODFT只有半数值频率,体系稍大一点就会极为吃力。对于第3步高精度的考虑标量相对论的单点,不是必须用相对论赝势,用高质量全电子基组并使用较好的标量相对论哈密顿(LUT-IOTC、DKH3等)也可以,只不过耗时会高不少而且精度也未必更好。