在赝势下做波函数分析的一些说明

 
经常有人问在赝势下能不能做AIM分析、Multiwfn (http://sobereva.com/multiwfn)支持不支持赝势等等问题。此文专门就此问题进行探讨。
 

1 全电子相对论计算

对于重原子，从第四周期开始（K~Kr）相对论效应就已经体现出来了，第四周期之后就必须考虑相对论效应了。对于使用全电子基组时，一般通过DKH等方法做标量相对论计算，比如在Gaussian中用DKH2方法时只需写上int=dkh2就可以了。

平时常用的基组的收缩方式并不是为了相对论计算而优化的，用于相对论计算并不适合，建议用相对论计算专用基组。比如(aug)-cc-pVnZ-DK系列比(aug)-cc-pVnZ系列明显更适合DKH2计算。有些人也将已有的非相对论基组重新调整收缩方式使之适合相对论计算。完全未收缩的基组也可以直接用在相对论计算上，但是计算量往往过大。

适合相对论计算的全电子基组有这些：

• (aug)-cc-pVnZ-DK系列：(aug)-cc-pVnZ的专用于DKH计算的重收缩版本。主要适合后HF计算，由于高角动量函数多，广义收缩程度又高，DFT计算用它的话亏。支持前六周期元素（K/Rb/Cs、Ca/Sr/Ba除外）和锕系的Th/Pa/U。
• def2系列重收缩版：def2系列基组从H到Xe都有全电子基组版本，但没有对相对论计算而优化。ORCA的一伙人对这系列基组进行了重收缩使之适合DKH2、ZORA计算。这也是ORCA对于H~Xe的DKH/ZORA计算时标配的基组。很适合DFT计算。基组定义可以令ORCA程序用PrintBasis关键词打印出来放到别的程序中使用。
• SARC（Segmented all-electron relativistically contracted）：ORCA程序的主要开发者提出的，是很适合DKH2、ZORA结合DFT计算的全电子基组，基组不太大结果也挺好。目前公开发表的只有镧系、锕系、第三周期过渡金属和Tl~Rn。很适合搭配def2系列重收缩版，这就可以涵盖周期表绝大多数元素了。后来又提出了SARC2，比SARC精度更高，但覆盖的元素目前还较少。
• ANO-RCC：元素覆盖全面，很适合考虑电子相关作用的相对论计算。但是基组巨大，收缩度特别高，高角动量函数又多，对于HF/DFT这样不需要很高角动量的方法来说就非常浪费了，稍大点的体系很难算得动。
• Turbomole的x2c-基组：def2的作者对标量和二分量X2C计算专门优化了全电子2-zeta和3-zeta基组，最高支持到Rn。包括x2c-SV(P)all、x2c-SVPall、x2c-TZVPall、x2c-TZVPPall，质量与def2系列同档次版本相仿佛，并且也定义了/J辅助基组。质量很理想
  
• WTBS：元素涵盖全面，是极小基，收缩程度很大，是给非相对论HF计算专用的，能接近HF能量极限。由于不是为相对论而优化的收缩的形式，用于相对论时收缩形式必须得改改才行。
• UGBS：元素涵盖广，由于完全没有收缩所以虽然不是特意给相对论计算提出的但也毫无问题，GTF虽不多但是基函数很多，计算费时，在Gaussian中还很容易出现积分精度不够的错误（需要用较高数目径向积分格点来解决，如int=grid=400434关键词）。此基组可以加上极化和弥散函数（见Gaussian手册）。
• Sapporo全电子相对论基组：札幌的一批人提出了DZ、TZ、QZ级别的用于非相对论和DKH3相对论计算的片段收缩基组，可在此获得：http://sapporo.center.ims.ac.jp/sapporo/Order.do。
• Hirao等人在JCP,115,4463(2001)提出了专供DKH3相对论计算的基组，参数通过最小化原子DKH3标量相对论计算得到的能量得到，覆盖1-103号元素，定义可在此获得：http://www.riken.jp/qcl/publications/dk3bs/periodic_table.html。
• Dyall搞了专门适合四分量、二分量计算的基组，是完全非收缩的，从第四及之后周期有定义，DZ、TZ、QZ级别都有。可在此获得：http://dirac.chem.sdu.dk/basisarchives/dyall/。
• RPF-4Z是一种性价比很高的4-zeta级别用于相对论计算的基组。RAGBS（Relativistic Adapted Gaussian Basis Sets）号称是相对论计算时没有变分塌陷问题的基组中最小的一种，可在此获得：http://basis-sets.iqsc.usp.br/basis-sets/。
• Jorge课题组开发了一套DZ到6Z级别的基组，涵盖了周期表大部分，有用于普通计算和专为DKH2计算优化的两类基组。可在此获得：http://qcgv.ufes.br/downloadbasis.html。
•  DGauss系列基组是专做DFT的DGauss程序御用的基组。定义了从H~Xe，对非相对论LSDA计算优化得到，包含不同档次。只有B~Ne公开发表了，其它的内置于DGauss和Gaussian程序中（在后者中叫DGDZVP、DGDZVP2、DGTZVP，DG=DGauss）。不适合直接用于相对论计算。
  

有兴趣者可参考《从H到Lr所有元素的全电子波函数文件》（http://sobereva.com/235），其中用全电子基组结合DKH2计算了H到Lr所有元素，输入输出文件都提供了。

2 大核赝势和小核赝势

量子化学中通常指的赝势(ECP)就是用有效势来表达内核电子的势场，这样就只需要计算外层的，也即与化学过程密切相关的电子就可以了，大大节约了计算量。使用了等效体现相对论效应的赝势(RECP)还避免了处理重元素时需要做标量相对论计算，量化研究中常用的赝势都是RECP。简单来说，如果被赝势代替的内核电子（赝化的电子）较少就叫小核赝势，用这样的赝势计算时，价层、次价层乃至更内层的电子都需要明确计算。如果赝化的电子数较多，就叫大核赝势，用这类赝势时只有价层电子是要明确地计算的。

显然，小核赝势比大核赝势耗时，原则上说要比大核赝势好，但是实际上也并非总如此，比如对于镧系锕系，记得Xiaoyan Cao就指出用DFT时小核赝势（及全电子基组时）结果还不如用大核赝势，可能是DFT都没对这么重的元素优化。如果能确定价态的话，DFT最好使用相应价态的大核赝势（有些赝势对于原子在不同价态时是不同的）。

对于一个不熟悉的赝势，想知道是大核赝势还是小核赝势，一种方法是看原文，另一种方法就是到BSE（https://www.basissetexchange.org），如果能查到相应赝势，从中查看有多少电子被赝化。比如对于SBKJC赝势的碘，用Gaussian格式查看赝势定义，就会看到这一行
I-ECP     3     46
这代表46个内核电子被赝化了，由于碘是53个电子，所以明确考虑的有7个电子，显然就是5s 5p价层上的7个电子，因此，这是大核赝势。而对于cc-pVDZ-PP，对于碘被赝化的是28个电子，因此比SBKJC时要多考虑18个电子，自然这就是4s 4p 4d上的2+8+10=18个电子（不考虑4f，因为其能量较5p更高）。对于碘，由于cc-pVDZ-PP不仅考虑了价层的电子还考虑了次价层的电子，所以是小核赝势（具体来说是小核赝势基组，见下文）。

Lanl系列赝势中，对于过渡金属，Lanl1赝势就是大核赝势，比如对于Fe，只考虑了4s 4d。Lanl2赝势就是也考虑了次价层电子的小核赝势，比如对于Fe，考虑了3s 3p 4s 4d。特别要注意的是，lanl1、lanl2对于主族元素（不包括第一和第二主族，下同）的定义是完全一样的，且是大核赝势！顺带一提，和多数赝势不同的是，Lanl赝势对第四周期元素并没有考虑相对论效应，Kr之后的才考虑了。

SBKJC赝势对主族是大核，对过渡金属是小核。

Stuttgart (SDD)赝势对多数过渡金属是小核，对一些过渡金属尤其是稀土金属元素也有大核版本。在Gaussian中使用时，可以用XYn形式来指定，n代表多少内核电子被赝化，可用的组合方式见手册Pseudo关键词末尾部分。Stuttgart赝势对于主族也有大核和小核赝势之分，不过Gaussian内置的只有大核的。

赝势和赝势基组是两个概念，只有二者结合才能用于实际计算。lanl系列、SBKJC系列、Stuttgart系列等各种赝势都有官方标配的赝势基组，也有一些第三方的赝势基组，比如(aug)-cc-pVnZ-PP实际上不是赝势的名字，而是基于Stuttgart赝势优化出的适合于相关计算的赝势基组。(aug)-cc-pVnZ-PP对于所有元素，无论过渡金属还是主族，结合的都是小核版本的Stuttgart赝势。def2系列基组，诸如def2-TZVP，对于前四周期都是全电子基组，但是从Rb开始，尽管还是叫做原先的名字，但本质上就成了赝势基组了，结合的赝势都是Stuttgart小核版本的赝势。

由于绝大多数赝势对主族元素都是大核赝势，如果想用小核赝势算主族元素，推荐用(aug)-cc-pVnZ-PP或def2系列，它对过渡金属和主族元素都是小核赝势。用Stuttgart小核赝势结合标配基组算主族也可以，但是Gaussian里没自带，需要去Stuttgart赝势的主页上拷贝下来：http://www.theochem.uni-stuttgart.de/pseudopotentials/clickpse.en.html。

3 赝势下能做波函数分析乎？

原理上说，只要严格满足了两个前提，赝势下做波函数分析就完全没问题：1 关心的只是价层或更外层的性质 2：只有很内层的电子被赝化，价层的电子基本不受影响。

举几个例子，假设前述条件完全满足：
1 分析静电势时，由于实际体系中原子内核的密度基本不受化学环境的影响而是球对称分布的，因此它对分子价层及以外区域的静电势的贡献会与等电量的核电荷产生的贡献精确抵消，所以把内核电子赝化不会影响分子价层及以外的静电势的分析。
2 分析键级时，由于内层轨道不参与成键，因此显然可以放心地把内层电子赝化。
3 计算原子间离域性指数时，由于内核电子定域性很强，几乎不会离域到其它原子上，因此赝化内层电子不影响结果。
4 通过等值面分析电子定域化函数(ELF)时，如果赝化内核电子，尽管会导致看不到内核电子的定域化域，但由于价层波函数仍保持实际情况，原子间共价键和孤对电子的定域化域依然照常存在。
5 分析AIM临界点、键径时，由于赝化了内核电子，使得原子核位置的电子密度的峰消失，会导致找不到核临界点(NCP)（就算找到也只是数值原因产生的巧合）。而且由于密度不再是简单地从核中心单调下降而变得略复杂，因此在内核附近区域还可能会产生一大堆额外的临界点混淆视听。但是对应成键的键临界点(BCP)，以及环和笼的临界点(RCP和CCP)的位置和性质并不会因此发生变化，完全可以照常分析。而寻找键径时自然不可能再找到连接核中心NCP和对应成键的BCP的键径，不过诸如BCP-RCP这样原本只在价层出现的键径则不会受到影响。

上述这些都是在满足那两个前提的理想条件下的结果。但实际情况更复杂，因为内核电子并不只分布在内核区域，它们还与价层区域电子存在重叠，这是因为用作基函数的高斯函数是平滑衰减的，内核区域的高斯函数因此会在价层区域露着“尾巴”。因此，对于基于实空间的波函数分析，比如分析密度、ELF等，赝化内层电子会多少影响到价层性质的定量乃至定性结果。即便是基于希尔伯特空间的分析，比如计算Mayer键级，一些价层以内的电子，尤其是次价层电子，也可能一定程度地表现出价层电子的性质，涉及到原子间相互作用，因此将它们赝化会影响分析结果。

由于次价层电子最接近价层电子，所以为了尽可能少地避免赝势对波函数分析结果的影响，次价层电子总应当保留，也就是总应当使用小核赝势，否则可能定性正确的结果也得不到。但是也不代表用小核赝势就保险了，赝化比次价层更内层的电子依然可能对价层的分析带来一定定量的扰动，但起码此时定性正确的结果可以获得。在JCC,18,416中，作者分析了赝势对AIM分析的影响，发现使用小核赝势时虽然价层的动能密度、拉普拉斯值与实际有些偏差，但结果还是可信的，该有的临界点也都有，而使用大核赝势时，有的BCP就不见了。

4 wfx文件与EDF信息

wfn文件格式是最早Bader小组开发的AIMPAC程序引入的，用于记录分子坐标、高斯函数和轨道波函数信息，已被Gaussian, GAMESS-US/UK, Firefly, Q-Chem等众多从头算软件所支持。这种格式的三大弊端是：
1 高斯函数最高角动量只支持到f。虽然其格式实际上也可以用于记录更高角动量，但是更高角动量的标号没有统一定义（不过从G09 B.01开始输出的wfn文件已经支持更高角动量了，Multiwfn也可以正确载入）
2 没有明确记录轨道自旋类型，是alpha还是beta轨道有时得靠猜
3 数据的记录精度有限

Gaussian从Gaussian09 B.01开始支持了wfx格式，用out=wfx即可输出，它是对wfn文件的扩展。wfx格式搞得很繁复，我个人认为格式定义得不好，程序读取起来也麻烦，但是带来了几大好处：
1 高斯函数角动量可以支持到无限高，标号有统一的定义
2 轨道自旋类型有了明确的记录
3 数据的记录精度得到了提高

除了在上述三点的改进，wfx格式最关键的是加入了一个EDF(electron density function)字段，这对于在使用赝势时进行基于密度的分析十分有用。这个字段记录了赝势计算时内核电子的分布，将这部分密度和实际计算得到的价层的密度叠加起来和全电子计算时的密度是几乎一致的。各元素内核部分的密度是内置在Gaussian程序中的，是事先通过多个s型高斯函数对全电子相对论获得的内核密度拟合得到的，会在Gaussian计算完毕后输出wfx文件时直接写入EDF字段。关于拟合有关的更详细信息见JPCA,115,12879。EDF字段定义的内核密度和赝势种类无关，是看元素和被赝化的电子数。比如对于Cu使用lanl2DZ和SDD，由于都是小核，所以.wfx文件里的EDF字段的信息是相同的。

注意：对于极个别体系，Gaussian某些版本生成的wfx文件的EDF字段信息不对。比如Ti原子在lanl2下本该有10个内核电子，但是在个别含Ti体系中，G09给出的wfx文件里对应的EDF字段只描述了4个电子，导致此错误的原因目前不明（在C.01版仍存在）。如果不放心，可以在使用Multiwfn分析wfx文件时，先进入主功能100里的功能4，然后选1，这会对整个空间的电子密度进行积分。如果积分的结果和总的电子数（内核电子+价电子）比较接近，那么说明内核电子的密度都在EDF字段中正确描述了。

Multiwfn支持wfx格式，并且会自动利用其中EDF字段的信息把内层电子密度在计算总密度（包括各种基于密度的函数，但静电势除外）的时考虑进去。不过从头算程序对wfx支持得还很不普遍。Firefly程序也支持产生.wfx文件，但是并不会把EDF字段写进去。

欲更详细了解wfn格式见《高斯fch文件与wfn波函数文件的介绍及转换方法》（http://sobereva.com/55）。欲了解wfx格式的具体定义参见http://aim.tkgristmill.com/wfxformat.html（可能需要fan wall）。

Multiwfn从3.4版开始有了一个极为重要改进，就是内置了Wenli Zou等人构建的EDF内核电子密度数据库。如果载入的.wfx文件里有EDF字段，那么Multiwfn照常从.wfx文件里读取EDF字段；而当.wfx文件里没有EDF字段，或者用的是其它包含波函数信息的输入文件（如.wfn、.fch、.molden等），默认情况下Multiwfn对于使用了赝势基组的原子会自动从自带的EDF数据库中选取合适的EDF信息载入，从而提供对内核电子的描述。而且这个Multiwfn内置的EDF数据库比Gaussian自带的EDF数据库的质量更好，支持的元素更多。因此，对于Multiwfn来说，使用wfx文件在描述内层密度上并没有任何优势。

如果不想让Multiwfn从.wfx文件中读取EDF字段，或者不想从自带的EDF库中自动读取EDF信息，可以用settings.ini文件里的readEDF和isupplyEDF参数修改规则，详见此文件中的注释。另外，Multiwfn中也允许从原子的.wfx文件中载入EDF信息。详见Multiwfn手册附录4，不过这种做法意义不大，因为靠内置的EDF库就足够了。

注意EDF信息提供的只是内层电子的密度，在Multiwfn中也只有直接依赖于电子密度的实空间函数，比如电子密度、电子密度梯度模、电子密度拉普拉斯值、约化密度梯度(RDG)等的分析会受到EDF的影响，而其它的实空间函数，比如ELF、自旋密度、静电势，以及不涉及实空间函数的波函数分析，如Mayer键级、Mulliken分析，都不会受到EDF的任何影响。

虽然EDF可以让赝势下的与电子密度有关的分析和全电子基组时符合较好，但是，笔者依然鼓励用小核赝势，因为在考虑EDF下，小核赝势的结果还是会比大核赝势更接近全电子基组的结果。

5 总结

综上所述，如果在赝势下做波函数分析，一定要用小核赝势。如果用了小核赝势，但还是还想获得更高精度结果，或者怀疑结果可能有误，那么不如干脆直接上全电子标量相对论计算（对第四周期元素，忽略相对论效应用普通基组计算也可以）。如果是在赝势下分析和密度相关的问题，如根据AIM理论对电子密度做拓扑分析、分析电子密度拉普拉斯值，用Multiwfn可以非常放心，因为无论哪种格式的波函数文件，程序都会自动读入EDF信息，从而把内核密度如实地描述出来，结果会和全电子基组时相符很好。