赝势的函数形式以及在量子化学程序中定义的方式

文/Sobereva @北京科音  2013-May-29

1 前言

赝势在量子化学计算中使用得极为普遍，有时需要用的赝势在量子化学程序中没有内置，就得通过自定义的方式添加，例如Stuttgart赝势对于Au按照Gaussian程序格式的写法是
AU     0
AU-ECP     4     60
g-ul potential
  1
2      1.000000000            0.000000000      
s-ul potential
  2
2     13.205100000          426.709840000      
2      6.602550000           35.938824000      
p-ul potential
  2
2     10.452020000          261.161023000      
2      5.226010000           26.626284000      
[...略]

虽然定义的格式对于各个量化程序都大同小异，大家也经常接触到，但是其含义很多人不很明白。一方面是主流的量化书籍中往往对赝势都不怎么详谈，另一方面网上也没有什么资料专门讨论这个问题。本文主要目的就是以不很严谨但是浅显易懂的语言说明白通常用的赝势的函数形式，特别是在量化程序中是怎么靠类似上述内容定义出来的。在此之前先谈一下赝势的基本特点。

2 赝势的基本特点

本文说的赝势是量化研究中最常用的赝势，如Hay和Wadt提出的Los Alamos赝势（或称Lanl赝势）、Stuttgart/Dresden赝势（或称SDD赝势）、SBKJC赝势（也称CEP赝势）等等。

赝势的主要用处在于：(1)将化学上不感兴趣的内层电子以等效的势场描述，从而不需要将内层电子显式地表达出来，大大节约了计算量 (2)可以等效地体现相对论效应。从第四周期相对论效应开始显现，但不考虑也问题不大；而对第五周期及之后相对论效应是不可忽视的，不考虑甚至结果定性错误。

使用赝势后，所得的轨道可以称为赝轨道。最低的一批赝轨道的能量和全电子计算时各个价层轨道能量接近或严格相一致；轨道形状上，在内层区域，这些赝轨道没有了实际价层轨道的节点，而在价层及靠外的区域（即大于cutoff半径rc的区域）赝轨道和实际轨道的形状是很接近或严格一致的。只要赝势设定得合适，靠赝轨道一般可以合理表现出实际价层轨道所展现的性质，如成键方式、价层区域电子密度分布等。

赝势一般是对于元素来定义的，有的时候还根据特定的价态、考虑相对论的方式、被赝化的电子数等因素来具体分为不同版本。对于分子体系，赝势就是体系中各个原子对应的赝势的加和。所以这也要求赝势要有较好可移植性，适用于不同化学环境。

赝势的定义是依赖于被作用的波函数的角动量的（即依赖于角量子数l）。而对于主量子数n和磁量子数m的依赖在一般的研究中是不考虑的。（2、4分量相对论计算时所用的赝势还依赖于总角量子数j，本文不涉及）

构建赝势的方法各有不同。例如Los Alamos赝势是形状一致型赝势的典型，其构建的方法是：(1)对指定的原子，选定好哪些轨道作为内核轨道，哪些作为价轨道，并且选定原子所处的组态 (2)对原子做数值Hartree-Fock全电子计算得到价层轨道波函数 (3)把价层轨道波函数的r<rc部分的节点抹掉而成为平滑函数（以多项式表达），r>rc区域的轨道波函数则不作改动。这样就构建出了赝轨道 (4)对于不同角动量的赝轨道，根据赝轨道波函数及其能量反推出一个势函数（暂以列表方式数值表达），使得赝轨道正好是HF方程外加这个势函数下的本征函数，且相应的本征值恰为实际价轨道的能量 (5)将数值描述的势函数以高斯函数拟合来表达。由于势函数是平滑的，所以这容易做到 (6)选取一定数量的原始高斯函数(GTF)，然后优化它们的指数和收缩系数，使得它们能尽可能好地表现赝轨道，这就定义出了与此赝势对应的赝势基组。
而以Stuttgart赝势为代表的能量一致型赝势的构建思路与形状一致型赝势有明显差别。它将目标数据选取为原子及其离子态的成百上千个电子态的总价电子能量（而非一个个价轨道能），通过不断调整赝势参数，使赝势与相对论全电子计算得到的这些态的总价电子能量的平均差异尽可能小。尽管此类赝势没有严格要求r>rc范围内赝轨道和实际轨道完全一致，但结果上看却也非常接近。

显然，用赝势就得用价电子基组。赝势在提出时都有标配的赝势基组。同一种赝势也可以有多种不同质量或不同研究者定义的赝势基组，如Los Alamos赝势有Lanl2DZ、Lanl2DZdp、Lanl08、LANL2TZ、LAV或LACV系列基组等等。对于Stuttgart赝势，除了它标配的基组外，(aug)-cc-pVnZ-PP系列、较重元素的def2-系列基组用的也都是Stuttgart赝势。如果没有经验，对于特定赝势，一般不建议使用未专门对它优化过的价层基组，不同赝势各自适用的赝势基组也不建议随意互换，这是因为不同赝势对应的赝轨道在内层区域径向形状是有差异的，乱用赝势基组会造成不能合理表现这部分区域。

3 赝势的函数形式与角动量问题

使用赝势时，体系的哈密顿为H = ∑[i]h'_i + ∑∑[i>j]1/r_ij，第一项是价电子的有效单电子哈密顿的加和，第二项是价电子间的静电相互作用。有效单电子哈密顿为h'=h+∑[a]U_a，其中h是常规的单电子哈密顿，即动能+核吸引势。U_a代表a原子对应的赝势。

有了H，就可以基于赝势求解化学问题了。赝势对于HF/DFT/post-HF计算都可以用。对于HF（或DFT）而言，HF方程只考虑价层电子，fock算符中包含只由价层电子构建的库仑项和交换项、∑[a]U_a以及核吸引势。这样求解HF方程得到的就是赝价层轨道。

下面来看看赝势的函数形式。由于赝势都是球形平均的，所以只需考虑径向坐标。各种赝势都写为这样的通式：
U(r)=U_L(r)+∑[U_l(r)-U_L(r)]*P_l  （加和是对于角量子数l，从0到L-1）
式中L往往取被赝化的内核轨道中的最高角动量+1。P_l=|l><l|代表l角动量投影算符，可以进一步写为对应的磁角动量投影算符的加和P_l=∑[m=-l→l]P_l,m。比如P_2投影算符投影出的是d成份（l=2），若将它作用在p轨道（l=1）上，结果就为0；若作用在d轨道上，则轨道保持原样；若作用在一个复杂的轨道上，则得到径向、角度方向可分离变量形式R(r)*∑c_lm*Y_lm(θ,φ)，其中R是此轨道的径向部分函数，Y_lm是l,m型球谐函数，c_lm则是此轨道中Y_lm球谐函数的分量。

上述赝势U(r)表达式表现的含义是对于任何轨道都施加U_L势。而对于轨道中的角量子数l在0到L-1之间的成份，则进一步施加相应的校正势U_l-U_L。换句话说，轨道中角动量l>=L的部分感受到的是相同的U_L势，而角动量l<L的部分感受到的则是依赖于l而因此各自不同U_l势。

前面已提到赝势是依赖于角动量的。假设内核轨道中角动量最高的为t，那么，对于价层轨道l=0,1,2...t角动量成份都必须单独定义相应的赝势。这是由于不同角动量的内层轨道与价轨道的相应角动量成份之间有着明显不同的相互作用（具体来说，是正交以及交换势）。对于l>t的情况，随着l增大，l与l+1的势的差异很快就变得十分微小，这是为什么在通常的赝势中，会选取一个L，令l>=L的角动量成份都使用相同的势U_L，而不再对它们单独定义势。L的选取使情况而定。例如，对于第一周期过渡金属，Hay和Wadt发现d和f轨道对应的势差异比较显著，而f和g的势就没什么差别，所以对于第一周期过渡金属Los Alamos赝势中对于f及以上角动量所用的赝势一致，即L=3，而对s、p、d都单独定义了赝势。而对于第三周期过渡金属则需要明确考虑更高角动量，Los Alamos赝势将L设到了4 (g)，Stuttgart赝势甚至设到了5 (h)。

注意赝势基组的最高角动量和赝势定义中的L没必要一致，一般L是大于或等于赝势基组中最高角动量+1的。比如对于Cu，Los Alamos赝势的L=3，而Lanl2DZ赝势基组中Cu的基函数的最高角动量是d，用以描述价层d电子。之所以基函数没有f但是赝势中却定义了角动量大于等于f成份所感受到的势，是因为实际计算的是分子体系而非原子体系，赝势的作用对象不是当前原子的基函数而是赝轨道，其它原子的基函数都会参与赝轨道的构成，因此赝轨道的形状是复杂的，投影算符对它们是可以投影出大于等于f的成份的，这部分也应感受到赝势，所以赝势不能只给s,p,d成份来定义。但是，假设计算的只是Cu这一个原子，那么此时的赝轨道，换句话说也就是没有内部节点的原子轨道，最高角动量只能到d，而投影不出f及以上的成份，因此此时只需要定义s、p、d受到的赝势就够了，f及以上成份感受的赝势即使定义了也没用。

为了实际求解时方便，数值描述的赝势一般用多项式与高斯函数的乘积来拟合来得到解析表达式。赝势一般的展开形式都是统一的：
r^2*U(r)=∑[k]d_k*r^n_k*exp(-ζ_k*r^2)，等价于U(r)=∑[k]d_k*r^(n_k -2)*exp(-ζ_k*r^2)。
其中r是径向坐标，d_k是展开系数，n_k可以为0、1或2，ζ是高斯函数的指数。正是由于赝势的实际表达形式如同高斯型基组一样遵循统一的规范，所以各种赝势原则上都可以在各种支持赝势的量化程序中使用。通常赝势表达方法是先对U_L(r)用上式进行拟合，得出对于l>=L角动量成份的赝势的d、n、ζ数据，然后对于l<L的角动量部分，分别按上式拟合出对U_L赝势的修正量U_l(r)-U_L(r)的d、n、ζ数据。

当量子化学程序使用某赝势时，必须有能力做到L角动量基函数的积分。例如GAMESS-US的积分代码支持的基函数最高只到g，而对于Hg，Stuttgart赝势的L=5 (h)，所以尽管其赝势基组的最高角动量只是d，在GAMESS-US中也没法用。注意输出的波函数中最高角动量函数与赝势的L无关。比如标准的wfn文件，最高只能记录到f角动量高斯函数，对于含Hg体系并且用了Stuttgart赝势的情况，由于赝势基组自身没有高于f角动量的基函数，若其它原子所用的基组也没有高于f角动量的基函数，那么wfn文件就可以正常记录此体系的波函数信息。对于波函数分析程序，诸如Multiwfn，支持的高斯函数最高角动量是h，用于这样的体系的分析也是完全没有问题的。

4 赝势的定义格式

赝势在不同文献中的表达方式、在不同量子化学程序中定义的格式都是大同小异的。下面按照Gaussian程序中定义赝势的格式来进行说明。赝势基组的定义方式和普通的基组定义方式没有任何区别，所以这里不提及。

下面这是Zn的Los Alamos赝势的定义，//后面是注释：
ZN     0       //ZN是被定义赝势的元素名，如果只对某个原子定义，则可以写原子序号。0是必须写的
ZN-ECP     3     18  //赝势名字；赝势的L值，比如此例对于f及以上角动量用相同的势，故L=3；用赝势描述的内核电子数，这里为18。赝轨道中不同角动量成份所受的势分别写在接下来的一个个block里。第一个block定义的是l>=L成份所受的势（即U_L），接下来依次定义l=0、l=1、l=2...直到l=L-1对应的势相对于U_L的修正量
f-ul potential   //从这里开始是U_L势的block。这一行是block的名字，可以随意设，不影响结果
  5  //U_L展开为前述的∑[k]d_k*r^(n_k -2)*exp(-ζ_k*r^2)形式，k=1、2、3、4、5
1    386.7379660            -18.0000000   //n_1、ζ_1和d_1
2     72.8587359           -124.3527403   //n_2、ζ_2和d_2
2     15.9066170            -30.6601822   //n_3、ζ_3和d_3，以此类推       
2      4.3502340            -10.6358989        
2      1.2842199             -0.7683623  
s-ul potential   //这个block定义的是s角动量成份所受的势相对于U_L的修正，即U_0-U_L。这行是block的名字，也可以随意设
  5  //U_0-U_L也是作如上形式的展开，展开为以下5项。下面每一行都是相应项的n、d和ζ。
0     19.0867858              3.0000000        
1      5.0231080             22.5234225        
2      1.2701744             48.4465942        
2      1.0671287            -44.5560119        
2      0.9264190             12.9983958      
p-ul potential   //和上一个block类似，此block定义的是p角动量部分对于U_L的修正，即U_1-U_L
  5
0     43.4927750              5.0000000        
1     20.8692669             20.7435589        
2     21.7118378             90.3027158        
2      6.3616915             74.6610316        
2      1.2291195              9.8894424        
d-ul potential   //和上一个block类似，此block定义的是d角动量部分对于U_L的修正，即U_2-U_L
  3
2     13.5851800             -4.8490359        
2      9.8373050              3.6913379        
2      0.8373113             -0.5037319       

由此例可见，每一个block的名字都是随便起的，比如d-ul也可以写为d and up，s-ul可以写成s-d，或者瞎写一个名字，诸如写写Mio_akiyama也无妨。关键的是，block的顺序不能错，一定是先定义U_L，再按角动量由低到高依次定义U_l-U_L，l=0,1,...L-1。

5 修改赝势的写法

前面已经说明了赝势的形式和定义方法，为了加深理解，这一节我们对赝势的写法进行各种修改来其考察影响。

下面是Pt的Stuttgart赝势
PT     0
PT-ECP     5     60
h-ul potential
  1
2      1.000000000            0.000000000      
s-ul potential
  2
2     13.428651000          579.223861000      
2      6.714326000           29.669491000      
p-ul potential
  2
2     10.365944000          280.860774000      
2      5.182972000           26.745382000      
d-ul potential
  2
2      7.600479000          120.396444000      
2      3.800240000           15.810921000      
f-ul potential
  1
2      3.309569000           24.314376000      
g-ul potential
  1
2      5.277289000          -24.218675000   
值得一提的是这个势对于大于等于h角动量的成份其实是没有势的作用的，因为其block只有一项，而且系数d为0，这一项其实只是走个形式。之所以对>=h角动量没有定义赝势，这和前面提到的Stuttgart赝势生成方式有关，因为定义>=h的赝势对计算原子各种组态的总价电子能量的影响是可忽略的，于是就不去弄相应的势了。

如果我们计算的是单个Pt原子，由于价轨道只由自身的s、p、d基函数来描述，没有f及以上角动量成份，所以赝势定义简化为下面这样，即去掉f、g以及>=h的定义，对Pt单原子计算结果是没有任何影响的：
PT     0
PT-ECP     3     60   //L=3，以下依次定义>=f的势以及s、p、d对势的修正
f-ul potential
  0   //含0项，>=f的成份将不受赝势作用。尽管计算Pt单原子考虑f是多余的，但按照基组定义格式，必须先定义>=L所受的势，所以这里走个形式
s-ul potential
  2
2     13.428651000          579.223861000      
2      6.714326000           29.669491000      
p-ul potential
  2
2     10.365944000          280.860774000      
2      5.182972000           26.745382000      
d-ul potential
  2
2      7.600479000          120.396444000      
2      3.800240000           15.810921000   

如果我们计算的是含Pt化合物，做将原始的赝势改为如上这样就会一定程度影响结果，因为赝轨道比原子轨道复杂得多，因此可以投影出>=f的成份。对于那些最高角动量有比较大限制的程序，比如GAMESS-US只能做到g，若非要用这样L达到h（或更高）的赝势，就不得不人为地修改赝势来对角动量定义进行截断（也就是让赝势定义只包含>=g,s,p,d,f的定义，且让>=g部分不受赝势作用），尽管计算分子时这会一定程度上影响结果。并且非要这么做的话需要在文章中注明用的是自行对角动量进行截断过的赝势。另一种办法就是换其它的赝势，比如对于Pt，Los Alamos赝势和SBKJC赝势的L=4，比Stuttgart赝势低一阶，所以GAMESS-US中就能够直接用。


再举一个例子，对于前面的Zn的Los Alamos赝势，我们进行改写，不设U_L势（即令之定义为空），且将原先U_L势和各个角动量的修正势U_l-U_L直接合并，得到下面的形式：
ZN     0
ZN-ECP     3     18
f-ul potential
  0       
s-ul potential
  10
1    386.7379660            -18.0000000        
2     72.8587359           -124.3527403        
2     15.9066170            -30.6601822        
2      4.3502340            -10.6358989        
2      1.2842199             -0.7683623
0     19.0867858              3.0000000        
1      5.0231080             22.5234225        
2      1.2701744             48.4465942        
2      1.0671287            -44.5560119        
2      0.9264190             12.9983958        
p-ul potential
  10
1    386.7379660            -18.0000000        
2     72.8587359           -124.3527403        
2     15.9066170            -30.6601822        
2      4.3502340            -10.6358989        
2      1.2842199             -0.7683623
0     43.4927750              5.0000000        
1     20.8692669             20.7435589        
2     21.7118378             90.3027158        
2      6.3616915             74.6610316        
2      1.2291195              9.8894424        
d-ul potential
  8
1    386.7379660            -18.0000000        
2     72.8587359           -124.3527403        
2     15.9066170            -30.6601822        
2      4.3502340            -10.6358989        
2      1.2842199             -0.7683623
2     13.5851800             -4.8490359        
2      9.8373050              3.6913379        
2      0.8373113             -0.5037319   
对于Zn单个原子的计算，这种赝势写法和原先的赝势的写法计算结果是完全一致的，因为s、p、d成份感受到的赝势与原先相同。但是，对于分子体系计算，由于赝轨道存在f及以上角动量成份，我们又没设这部分的赝势，而原先定义中这部分赝势是非0的，所以计算结果将会有差异。