使用Gaussian做镧系金属配合物的量子化学计算

Using Gaussian to perform quantum chemistry calculations for lanthanide metal complexes

文/Sobereva@北京科音  2021-Jan-13

过渡金属配合物体系的量子化学计算比主族体系的计算通常更复杂，因为在判断自旋多重度方面需要更多考虑，SCF收敛整体更难（d轨道近简并所致），而且还有更高几率收敛到不稳定波函数。而镧系锕系配合物的计算甚至更困难，由于f轨道近简并导致SCF比过渡金属配合物还更难收敛（除非镧系锕系原子用大核赝势，此时不需要描述f电子），检测波函数稳定性的必要性更高，而且中心金属配位数也往往更高、考虑更复杂。本文给大家一个用Gaussian做DFT计算镧系金属配合物12-冠-4-五水合钬(III)阳离子的例子，希望通过此例令读者了解到镧系锕系金属配合物计算的相关知识和需要注意的相关问题，从而能举一反三地研究其它此类体系。

注：如果你的量子化学研究经验有限，本文提到的所有博文都要读，否则本文里一些说明难以充分、正确理解。如果你连用Gaussian算有机类体系都还不能熟练地计算，不建议现阶段看此文然后就开始一通折腾镧系锕系计算，应当看的是《谈谈学量子化学如何正确地入门》（http://sobereva.com/355）。

本文用到的计算程序是Gaussian 16 A.03，本文所谈的与Gaussian有关的情况对于这个版本都是适用的。相关分析用的是2021-Jan-10更新的Multiwfn 3.8(dev)版。

本文只涉及最常用的DFT计算。过渡金属配合物的泛函测试文章很多，但算镧系锕系用什么泛函最合适没有特别系统全面的测试。PBE0是值得优先考虑的，在不少文章里用于算镧系配合物，比如Inorg. Chem., 58, 411 (2019)里用就PBE0结合SDD赝势做了含Pr体系的计算。在SARC全电子基组的原文J. Chem. Theory Comput., 5, 2229 (2009)和Jorge 2-zeta全电子基组原文Chem. Phys. Lett., 643, 84 (2016)里，作者基于PBE0做全电子计算得到了挺好的结果。

PBE0不是唯一选择，笔者在《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）中提到的那些可以算过渡金属配合物的泛函也都可以尝试用于镧系配合物的计算，诸如BP86、TPSS、TPSSh等。一般体系计算非常常用的B3LYP不是上佳的选择，虽然在不少文章中也有人拿它算镧系。

镧系锕系元素相对论效应很显著，必须考虑相对论效应。其配合物可以用全电子基组结合相对论哈密顿（常见的有DKH2、DKH4、ZORA、X2C）做全电子相对论计算，也可以对镧系锕系原子用考虑了相对论效应拟合的赝势计算。在Gaussian中，对于这类体系的优化和振动分析一般应当基于赝势计算，因为Gaussian没有相对论哈密顿的解析导数，故这类任务极其耗时。

对赝势缺乏了解的话，先看《谈谈赝势基组的选用》（http://sobereva.com/373）和《赝势的函数形式以及在量子化学程序中定义的方式》（http://sobereva.com/188）。另外，笔者在北京科音中级量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）里有十分全面细致的关于赝势的专题讲解。对于计算镧系锕系，首选Stuggart赝势（亦称SDD赝势）。Stuttgart赝势对镧系锕系有大核和小核版本，在《详解Gaussian中混合基组、自定义基组和赝势基组的输入》（http://sobereva.com/60）中我专门介绍过。对于下文要计算的原子序号为67的Ho来说，大家去看Gaussian手册的Pseudo关键词部分（http://sobereva.com/g09/k_pseudo.htm），可知Gaussian对Ho内置了MWB28、MWB56、MWB57、MHF56、MHF57五种Stuttgart赝势。其中MHF的不用管，因为拟合赝势时没考虑相对论效应，结果根本没法用，诸如算的配位键键长都明显不合理。MWB28是以Wood-Boring方式考虑准相对论效应的相对论赝势，是小核赝势，把1s 2s 2p 3s 3p 3d共28个内核电子用赝势代替了。MWB56和MWB57是大核赝势，把1s 2s 2p 3s 3p 3d 4s 4p 4d 4f电子都赝化了。其中MWB56是专门给Ho(III)用的，Ho(III)相当于Ho丢了全部两个6s电子和11个4f电子当中的一个，由于1s 2s 2p 3s 3p 3d 4s 4p 4d一共有46个电子，再加上剩下的10个4f电子，因此一共是56个电子被赝化了。MWB57是专门给Ho(II)用的，其此时的全部11个f电子都被赝化了。用大核赝势的好处在于，由于在配合物中近简并的4f轨道不需要被描述，因此SCF收敛明显比较容易，而且相对于小核赝势来说计算也更快。大核赝势的缺点是精度一般不如小核赝势，而且可移植性差，不适合研究配位情况显著改变的化学过程（如涉及到配位键形成和断裂的反应、异构化），还需要自己恰当判断原子的价态。小核赝势用起来明显更为稳妥、普适，但缺点是耗时更高，而且由于4f电子的近简并导致SCF极其难收敛（实际算过的人肯定都深有体会），这令耗时更高。我建议小核赝势用在几何优化结束后的单点能计算、各种属性的研究（如计算激发态）、产生用于做波函数分析的波函数等场合（见《在赝势下做波函数分析的一些说明》http://sobereva.com/156），以及用于配位情况明显改变的过程的研究。

赝势和赝势基组是两码事，Stuttgart赝势有不同的与之兼容的描述价电子的赝势基组。对于DFT计算且要求不太高时，就用Gaussian里自带的Stuttgart赝势标配的赝势基组即可。如果想让结果更好，对镧系可以用def2-TZVP或更高级的def2-QZVP，它定义了赝势基组部分，对应的赝势是Stuttgart的MWB28小核赝势。在撰文时def2系列对于锕系没有定义，但可以用对锕系有定义的def系列，它也是定义了赝势基组部分并结合Stuttgart小核赝势。注意Gaussian至少截止到撰文时最新的G16 C.01仍没有内置镧系锕系的def/def2基组，必须去Turbomole基组库当中拷，见《在线基组和赝势数据库一览》（http://sobereva.com/309）。Lanl2赝势+标配基组支持La、U、Np、Pu，但用得很少，没必要考虑。CRENBL大核赝势也支持镧系锕系，但如今鲜有人用。

Gaussian支持DKH2和DKH4标量相对论哈密顿计算（另外还支持的DKH0和RESC太烂，这里就不提了），DKH2非常常用，比起非相对论计算耗时仅高一丁点。DKH4在原理上更精确，耗时比DKH2也就多个几分之一。对这两种情况的基组建议： • 用DKH2的情况：我建议对镧系元素用Neese搞的SARC基组，这是很好的3-zeta档次的基组，不特别昂贵结果又挺好。SARC有针对DKH2和ZORA优化的两种版本，显然当前情况要用前者。想要对镧系结果更好还可以用SARC2的DKH2版，是4-zeta档次的。二者在BSE上都有，不知道BSE是啥的话看《在线基组和赝势数据库一览》（http://sobereva.com/309）。SARC对前四周期没有定义，因此对于配体部分需要用其它基组，非常适合的是用Neese他们搞的针对DKH2重收缩的def2系列基组，但这没有公开发表，而只是内置于他们的ORCA程序中，在ORCA程序中用PrintBasis关键词打印出来定义再整理成Gaussian的格式虽然可以用，但稍麻烦，BSE上也没法直接获取。另一个对配体可以考虑的是Jorge专为DKH2优化的一套基组，从DZ到6Z档次都有，DZ和TZ版覆盖1~103号元素，BSE上搜Jorge即可下载到定义。如果你是做后HF计算，配体可以用带-DK后缀的给DKH2优化的Dunning相关一致性基组。 • 用DKH4的情况：我建议用Weigend等人（def2系列基组的作者）提出的x2c-系列基组，涵盖前六个周期所有元素（包括镧系），在BSE上可以下载到基组定义，但在Turbomole基组库中种类最全。此基组专门给X2C相对论哈密顿计算而设计，用于精度也很高的DKH4计算也完全可以，此基组原文测试发现结合X2C还是DKH4误差差不多。x2c-有用于标量相对论计算（后缀为all）和二分量相对论计算（后缀为all-2c）的两类版本，后者更大，对于DKH4标量相对论计算用前者就够了。x2c-和def2系列一样有不同尺寸的版本，比如x2c-SVPall、x2c-TZVPall、x2c-QZVPall等。

相对论量子化学的原理和全电子相对论计算的各种相关知识我在北京科音（http://www.keinsci.com）的高级量子化学培训班里会很详细展开讲，并全面盘点全电子相对论计算可以用的所有基组，本文就不再多说了。

下面就演示具体计算例子12-冠-4-五水合钬(III)阳离子，涉及到的所有文件都可以在http://sobereva.com/attach/581/file.rar下载。

作为例子要计算的12-冠-4-五水合钬(III)阳离子在剑桥结构数据库(CSD)里面有，不用我们自己画。它对应GINREA条目的阳离子部分，除此之外还有氯阴离子、游离的水分子，本文我们不考虑后两者，计算在真空下当做孤立体系来进行。

进入CSD搜索页面https://www.ccdc.cam.ac.uk/structures/，Identifier(s)里输入CSD识别码GINREA，就可以进入下面的页面

我们要导出其坐标。没必要先下载cif文件然后再提取分子，直接导出网页中显示的结构即可。做法是在3D窗口中点右键，选下面的选项

之后会出现一个文本框，把里面的内容复制到一个文本文件里，改名为GINREA.mol。

用GaussView打开这个mol文件，看到下图。由于GINREA条目是是X光衍射测定的结构，氢的位置一般没法准确确定或完全不能确定，水上只有氧原子。相关知识见《实验测定分子结构的方法以及将实验结构用于量子化学计算需要注意的问题》（http://sobereva.com/569）。

把三个氯离子，以及没有配位的两个水的氧原子都删掉，然后点击GaussView的加氢按钮，给配位水的氧上面都点上两个氢原子。如下所示

将当前结构保存为optfreq.gjf。

这个配合物的自旋多重度取决于Ho(III)的10个4f电子的排布，可以为单重态、三重态、五重态。在几何优化与振动分析过程中，按照前面的建议，我们用大核赝势来做这一步，此时4f电子是被赝化的，所以Ho(III)上的价电子就只有将5s和5p轨道都占满的8个电子了，而配体部分都是闭壳层电子结构，因此计算时应当设成单重态来算。

将optfreq改写为如下内容。

%chk=optfreq.chk #p pbe1pbe/genecp em=gd3bj int=fine opt freq