电子结构理论与计算.pdf

电子结构理论与计算 李震宇 中国科学技术大学合肥微尺度物质科学国家实验室 Version:20101228 1主要内容 主要内容 第一章第一章 电子结构理论基础电子结构理论基础 BO 近似，HF 理论，电子关联与后 HF 理论，密度泛函理论（DFT），交换关联泛函，赝势，半经验电子结构理论 第二章第二章 电子结构计算方法电子结构计算方法 原子基组，平面波基组与周期性边界条件，自洽场（SCF）第三章第三章 电子结构分析电子结构分析 电荷密度，分子轨道，能量 第四章第四章 电子结构理论与计算的应用电子结构理论与计算的应用 动力学与统计力学性质 附录附录 量子力学基本原理，常见数值算法 2前言 前言 理念：1.实用：介绍计算参数的意义和选择策略；2.简洁：注重物理概念，不过分纠缠于细节。参考资料：1 Richard M Martin,Electronic Structure:Basic Theory and Practical Methods.Cambridge University Press,Cambridge,2004.2 Christopher J Cramer,Essentials of Computational Chemistry:Theories and Models.2nd Ed.Wiley,England,2004.感谢所有帮助我将本讲义电子化的同学，以及所有其他选修本课程的同学。讲义中错误和疏漏之处在所难免，如有任何意见与建议，请发送到，十分感谢！3第一章 电子结构理论基础 第一章 电子结构理论基础 由于计算能力较弱，早期的科学研究以理论和实验的互动为主，其中理论研究通常使用唯象理论。随着科技的发展，现在已经有可能完全从量子力学基本原理出发来计算体系的各种性质。这就使得计算成为一种重要的独立的科学研究手段，与理论和实验形成三足鼎立之势。图示：新的科学研究范式三者之间有紧密的联系：计算可以帮助理解实验结果，例如扫描隧道显微镜(STM)观测到的是局域态密度，要从中获得几何结构的信息，需要计算模拟来配合；计算也可以验证理论猜想，例如通过计算机实验来检验材料设计原理。我们希望通过科学计算获得理论设计与计算表征的能力。体系的大多数性质都决定于它的电子结构，即电子的空间与能量分布。1.1 基于波函数的量子化学方法基于波函数的量子化学方法 1.1.1 波恩波恩-奥本海默奥本海默(BO)近似近似 在原子中，原子核的质量比电子的质量大很多，运动也比电子慢得多。1 amu(12C/12)是电子静止质量的 1822.83 倍。因此，通常在研究某一瞬间电子结构时可以忽略原子核的速度，假定它和原子核长期固定在某位置时的电子结构一样。这就是 BO 近似。其核心思想是把核的运动和电子的运动分开处理：处理电子运动时，认为核是固定不动的；处理核运动时，认为快速运动的电子建立一个平均化了的负电荷分布，核在电子的负电荷平均场中运动。对于包含许多电子和原子核的体系，其定态薛定谔方程可写为：,HEr Rr R其中 222222222IJIIiIiiIijIJIiiIijIJe Z Ze ZeHmmrrr-引入 BO 近似，对原子核和电子进行分离变量,;Nelr RRr R elelelHERrRr其中 NR是描写原子核状态的波函数，它只和所有核的位置有关；而 elHR和4;elr R是描写电子状态的哈密顿量和波函数，以原子核位置为参数【输入参数：POSCAR】。E R和原子核的位置有关，通常被称为势能面(PES)。BO 近似又称为“绝热近似”，它认为电子总是不断地适应核位置的变化，或者说跟随着核且始终围绕着核来运动，运动过程中不会发生从一个势能面 E1至另一个势能面 E2的跃迁，即吸收或放出能量发生电子能级跃迁。图示：石墨烯能量随 C-C 键长变化的势能面1.1.2 单电子近似单电子近似 由于电子之间存在相互作用，所以当电子数较多时，即使在 BO 近似下薛定谔方程仍然很难求解，需要作进一步的简化。最粗糙的近似是假设电子之间没有任何相互作用。这样，体系的哈密顿算符简化为：iiHh其中 222IiiIiiIe Zhmr 于是，我们可以先只考虑体系中的某一个电子，其定态薛定谔方程为iiiih，本征函数i为单电子波函数，通常也称为分子轨道(MO)，本征值i为对应的分子能级。下面的问题是怎么从单电子波函数写出总的多体波函数呢？对无相互作用电子，考虑泡利不相容原理，体系的波函数及其对应的能量本征值为12.HPn5iiE这种由分子轨道直接相乘得到的多体波函数称为 Hartree 积(HP)。显然，它不满足电子（费米子）全同性所要求的波函数的反对称性。为简单起见，先讨论无相互作用两电子体系12Hhh其中 221112hvm r 222222hvm r假设 ,12ababE 是两电子体系本征方程的解，则 12122112abababP 也是体系本征方程的解。这样，对应于同一个本征值，有两个线性独立的本征函数，但这两个本征函数都不满足反对称条件，为此作一线性组合，得：112212Aabab为体系满足波函数反对称条件的解，用行列式的形式表示为 111222abSDab该行列式被称为 Slater 行列式(SD)。对于由 N 个电子组成的全同粒子体系，可以将上述结果进行直接推广，相应的定态波函数为 12121212121112221,!NNSDNNNNNNNr rr上式的反对称性是很明显的，因为任意两个粒子坐标的互换，相当于行列式中相应的两列元素的互换，将改变行列式的符号。同时，若有两个或两个以上单粒子波函数相等，则0SD，意味着这样的状态是不存在的，这就是 Pauli 不相容原理。1.1.3 平均场近似与平均场近似与 Hatree-Fock 自洽场自洽场 单电子近似并不意味着需要完全忽略电子间的相互作用。比如，在考虑单电子问题时，我们可以将体系中其他电子的作用近似看成一个平均场，并假设该场不随时间而改变。平均场近似下的单电子哈密顿算符可以通过对能量(elH在多体波函数下的期望值)变分得到。对6HP 多体波函数，我们有2222IiiiIiiIe Zhvmr 其中jijj iijvdrr2jje因为在求平均场时其他电子的坐标都是积分变量，所以ih仅与第 i 个电子的坐标有关。对这种采用 HP 多体波函数的平均场近似，体系的总能量是222|12ijHPiijiijijeEd drr r其中第二项为重复计数修正。我们看到，在平均场近似中，需要求解单电子本征方程才能得到体系的解，但要求解单电子方程又必须知道各个轨道的电荷分布。这看起来自相矛盾。在实际操作中我们可以从 N个初始单粒子波函数 0i出发，建立相应的(0)ih，然后通过解单电子本征方程，得到一组新的分子轨道 1i，再用这组 1i来建立新的(1)ih，重复这个过程，直到最后得到的 ni与上一步的1ni相等或近似相等。这样一个过程称为自洽场(SCF)方法【输入参数：EDIFF】。Fock 建议将 Hartree-SCF 过程推广到采用 Slater 行列式作为多体波函数的情形。对两电子体系的电子库仑相互作用算符，我们有【作业】12121SDSDababd dJKrr r其中 22121212121112abHPHPabJd dd drrr rr r 121211122abababKd drr r分别代表电子间的库仑排斥能和交换能。交换能 Kab是采用 Slater 行列式以后比 HP 多出来的同自旋电子间的一种相互作用。Roothann 提出可以基组表象来建立了 HFSCF 的矩阵代数。如果将分子轨道写成一些基函数的线性组合，则单电子哈密顿量在这组基组下可以表示为一个矩阵。常用的基组有原子轨道和平面波。对闭壳层体系，最后可以得到久期方程det(F-ES)=0，其中单电子哈密顿量f通常被称为 Fock 算符，其矩阵元为7221122IIeIFZPmr 其中 212121122ed dr r r2occiiiPa a由于 Fock 矩阵的表达式中包含密度矩阵，我们需要自洽求解。至此，我们通过如下步骤得到了作为量子化学基础的 HartreeFork 理论。1.1.4 电子关联与电子关联与 post-HF 方法方法 HF 理论采用单行列式多体波函数，它假定一个电子在由原子核和其它电子形成的平均势场中独立运动，没有考虑电子间的瞬时相关作用，也就是电子关联。但是，HF 理论中采用的 Slater 行列式波函数满足 Pauli 不相容原理，所以它在一定程度上包含了同自旋电子的关联作用，这部分电子关联被单独命名为交换相互作用。所以，通常所说电子关联能是指精确的基态能量与 HF 能量之差：Ecorr=E-EHF。关联能一般只占体系电子总能量的 0.32%，但是它可以与化学反应的反应热或者活化能相当甚至更大。一般说来，化学反应途径上电子关联程度不同。所以，在考虑化学反应和电子激发等问题时电子关联十分重要。由电子瞬时作用产生的电子关联称为动态关联。此外，还有非动态关联，它指体系几个对称性相同的电子组态函数（单行列式）接近简并时导致HF 方法给出很大的误差。常用的电子关联计算方法【输入参数：计算方法】有组态相互作用（CI）：ConfigurationInteraction耦合簇方法（CC）：Coupledclustermethod多体微扰理论（MPx）：MollerPlessetperturbationtheory密度泛函理论（DFT）：Densityfunctiontheory前面三种都是以 HF 的为基础，经常被称为 postHF 方法。而 DFT 基于完全不同的思路，将在下一节中专门介绍。1.CI 理论 可以通过引入激发组态来考虑电子关联800|.vir occvir occaaababiiijijaia b ijCCC 一般用 S,D,T,Q 等字母来表示从低到高的激发组态，其中最简单的情况是只引入单激发组态(CIS)。CIS 可用来计算激发能，但对基态能量没有修正【作业】。最复杂的情况是考虑所有激发组态，即 FullCI，它是特定基组下的精度极限，一般只能处理非常小的体系。考虑到精度与计算量的平衡，一种有效的方法是将分子轨道分为不同的类型。对最重要的分子轨道，如 Fermi 能级附近的轨道，考虑所有可能的电子组态。这些轨道组成空间称为全活性空间(CAS)。对次重要的轨道，考虑一定限制条件下的可能的电子组态(RAS)。在确定需要考虑的电子组态以后，通常在优化 CI 组态系数的同时优化单粒子轨道，即进行有多个组态的自洽场(MCSCF)。在 MCSCF 的基础上还可以再做 CI。这时，CI 不以 HF 单行列式为基础。而是由 MCSCF计算的波函数中选择组合系数较大的组态作为参考组态，从这些参考组态的占据轨道上移动电子到虚轨道上去，产生新的激发组态函数。由于有多个参考组态，这种方法称为 MRCI。2.Coupled-Cluster 方法 耦合簇理论是目前最流行的包含电子关联效应的量子化学方法之一。耦合簇理论通过对波函数施加一个如下的 ansatz 来改进对体系的解：TccHFe其中T为簇算符。如果将T写成123.TTTT其中2occvirababHFijijij a bTt，则 FullCI 可以表示为:(1)HFT。CCD：23222(1)2!3!CCDHFTTT9对有限大小的 Fermion 体系,簇算符泰勒级数的长度是有限的。由于簇算符是取幂的，簇算符中的一个单激发算符同时也激发其它的高阶态。3.MP 微扰方法 CI 理论基于量子力学的变分原理，将试探波函数写成不同电子组态的线性组合。量子力学中另外一种常用的方法是微扰论。Muller 和 Plesset 提出了一种以 HF 波函数为零阶波函数的微扰方法，称为 MP 微扰理论。多体微扰哈密顿量为的一般形式为(0)HHV其中，零阶哈密顿量为单电子 Fock 算符之和(0)iiHf(0)HF2,1()2occoccijijj ii jijeVJKr可以证明，HF 能量就是一阶 MP 微扰修正后的能量【作业】。同 CI 方法比，用小基组进行MP 方法计算并没有什么实用价值。经验表明，选择大基组低级别的 MP 计算结果比选择小基组高级别的 MP 方法好。最常用的 MP 方法为 MP2。与 MRCI 类似，我们同样可以在 MCSCF 的基础上再做微扰修正。最常用的是 CASPT2，其零阶波函数即为 CASSCF 得到的波函数。对上面提到的 PostHF 量子化学方法，一个计算精度的粗略排序为：HFMP2MP3CCDCISDCCSDMP4CCSD(T)。另外，高精度量子化学计算还有 G2、G3 序列，其思想为借助误差抵消提高精度。【调研：PostHF 的最新进展】1.2 密度泛函理论（DFT）密度泛函理论（DFT）传统量子力学范式：外势 v(r)多体波函数(.)rrr可观测的物理量(observables)。电荷密度也是一个可观测物理量，可以写为：*2322().(.)(.)nnnnN dddrrrrrrrrrr对单体算符我们有*1 21 212*1 211 212|()()()()()NiNNiNNNood ddNod ddd on rrrrrrr rrrrrrrrr rrrr1.2.1HohenbergKohn 定理定理1：多体理论：多体理论定理一：全同费米子系统非简并基态的密度 n 唯一地决定了外势。这样我们得到新的范式：n(r)v(r)observables证明：假设两个外势av、bv，两者确定同一个电子几率密度 n0，则：1 P.HohenbergandW.Kohn,Phys.Rev.1964,136,B8641000000000|ababbabbbbabbbEHHHHvvE 由于av和bv为单粒子算符，所以00()oaabbEd vv nErr同样有00()obbaaEd vv nErr两不等式相加：0000ababEEEE故外势、电子几率密度之间存在一一对应关系，即 n 可作为体系的基本量。定理二：给定外势 v，存在 Fn定义在所有非简并基态密度 n 上，使得能量泛函()()vE nd vnF nrrr当0nn，即取基态电子密度时，能量取得唯一的最小值。证明：（Levi 方法，分两步走）min|min min|minmin|()()min()()nnnnnEHHTUd vnF nd vn rrrrrr1.2.2KohnSham 方程方程2：有效单体理论：有效单体理论Levi 泛函可写成动能和势能两部分：F nT nU n动能的主要部分可以通过参考具有相同密度的假想的无相互作用体系得到：2*2()()2nSiiiT ndm rrrUn的主要部分为：1()()2|HnnUndd rrrrrr()()XCtotSHSHEETVUTTUU如果将电子几率密度表示为：2()|()|occiindrrr2 W.KohnandL.J.Sham,Phys.Rev.1965,140,A1133.11通过变分法可得到如下单电子方程：【作业】2()()()2extHXCiiivvvm rrr需要注意 DFT 引入的“单电子”近似和 HF 理论中的区别。在 DFT 中，KS 轨道的引入仅需要保证上面的电子密度表达式，并没有对应的多体波函数。在 HF 理论中，假定了 Slater 行列式形式的多体波函数，并在此基础上写出近似的总能，再做变分；而 DFT 中用来变分的系统总能原则上是严格的，所以，如果知道交换关联泛函的严格形式，DFT 也是严格的。这同时也意味着 KS 轨道并不是严格意义上的单电子轨道。虽然实际上 KS 轨道经常被用来作为近似的分子轨道，但是严格说来它除了计算电子几率密度外它没有其它物理意义，KS 轨道能量也不同于分子轨道能量。与 HF 比，DFT 对系统能量的计算精度大大提高，而计算量略有下降。1.2.3 交换关联泛函【交换关联泛函【输入参数】电子之间的交换项来自于 Pauli 不相容原理|SDSDXnnHEnUUn|SDSDnnSHXTUT nUnEn minmin ()|CSHXSDSDnnnnEnF nT nUnEnTUTU 显然，EC，所需的积分数目由约 N(N+1)/2 个|slijG vG变为 N 个|KBilmG。KB 赝势最大的缺点是它的形式不再保证微分方程式解里“朗氏根(Wronskian)”的性质，对于一个给定的量子数 l，有可能存在一个有节点的波函数比零节点的波函数能量还低。因而，这就允许了 ghost state 的产生。所幸的是可以通过一个简单的判据判断 ghoststate 是否存在PRB1990,41,12264，并通过参数的调控来加以避免的。不管怎样，使用前对赝势的仔细测试是非常重要的。1.3.3.超软赝势（超软赝势（USPP）与投影缀加波（）与投影缀加波（PAW）方法）方法对任意参考能量，非局域赝势可以写成|iiNLpsiiv，其中辅助局域波函数|()|psiilociTv。推广到多个参考能量，可定义|ijijB,1|()|ijijjB,|NLijijijvB 容易验证 ()|pspslocNLiiiTvv 且当模守恒条件|0pspsijijRijRQ 成立时，Bij和 vNL是厄米的。图：氧原子 2p 径向波函数（实线）、模守恒赝波函数（点线）、超软赝波函数（划线）。摘自Phys.Rev.B 1990,41,7892 由于模守恒条件条件的限制，NCPP 对第一周期元素和过渡金属元素并不能显著减低计算量。于是，人们提出了不遵循模守恒条件的超软赝势7。这时，需要定义 1|ijijijSQ ()|NLijjijijijvBQ 7 D.Vanderbilt,Phys.Rev.B1990,41,789216这样，我们有|pspsijRijRS ，且()|pspslocNLiiiTvvS 在 USPP 中，最后的电荷密度表达式与 Q 以及|有关。1994 年 Blchl 提出 PAW(Projector augmented wave)方法8，将全电子波函数和赝波函数通过以下变换联系起来：,psRRTTIT与 USPP 类似，PAW 方法也引入投影算子。在原子周围作分波展开：|pspsmmmmpsmmpsmmmmccTI 有了变换关系 T，通过求得赝波函数就可以获得体系的能量,电荷密度等信息。AT AT 1.4 半经验电子结构计算方法半经验电子结构计算方法半经验电子结构计算的主要思想是在就单电子哈密顿量时引入一些近似来提高计算速度。这里我们主要关注紧束缚（tightbinding,TB）近似。Slater 和 Koster（SK）提出，TB 近似可以作为一种插值方法，首先通过 Brillouin 区高对称点的计算结果来获得 Hamiltonian 中的参数，再通过这些参数对整个 Brillouin 区进行计算。为了将 Hamiltonian 参数化，他们做了如下简化：(1)正交化基组；(2)只考虑近邻或者稍远的原子之间的相互作用；(3)只考虑价电子；(4)在求 Bloch 和之间的 Hamiltonian 矩阵元时，忽略*3()()nimjHd rrRrR中的双中心以上的项。双中心积分可以通过方向余弦由一组 SK 积分求得Phys.Rev.1954,94,1498。如果只考虑 s,p,d 轨道，SK 积分的个数为 14 个：ss,sp,pp,pp,sd,pd,pd,dd,dd,dd,ps,ds,dp,dp。如果 SK 积分参数选择合适，TB 方法可以有较好的可移植性，使得研究不在参数拟合结构库中体系成为可能。1.4.1NRLTBJ.Phys.:Condens.Matter2003,15,R413在原始的 SK 模型中，基组被假设是正交的。这可以通过对原子基组做 Lowdin 正交化实现。如果不做此限制，则重迭矩阵不是单位矩阵，原来的矩阵本征值问题变成一个广义本征值问题。重迭矩阵的矩阵元与 Hamiltonian 矩阵元一样，也可以用 SK 积分来表示。引入重迭矩阵产生了更多的参数，使拟合更方便。同时 Lowdin 正交化过程通常使得基组更离域，从而最近邻近似的有效性变差，这也鼓励我们引入重迭矩阵。因此，在 NRLTB 方法中采用非正交基组。根据 DFT，系统总能量可以写成带结构能量和一个密度泛函的和33()()(2)nnd kEF nkr。8 P.E.Blochl,Phys.Rev.B1994,50,17953.17如果能带能量原点取在 Fermi 能级处，F 是排斥能。而当能带能量原点取在价带底，F 是吸引能。通过对能带进行一个平移()()()nnF nNerkk，可以将总能简化为带结构能量。偏移量与密度有关意味着 onsite 能量受环境影响。为此，对单元素体系可以定义2()IJRIIJJeF R，其中5()()1CCR RlRRF Re。这样 Hamiltonian 的 onsite 项可以定义为24233Illlllhabcd。对多分量体系，可以对每种不同的元素分开统计。SK 参数可以写为22()()()llhRllllllllhRefRgR eF R，重迭矩阵元可以写成类似的形式。这样，对单元素体系，共需拟合 93 个参数。为了保证重迭矩阵的正定性，一个改进的重迭矩阵元表达式为223()()()llsRllllllllllsRpRqRrR eF R1.4.2DFTBPhys.Stat.Sol.B2000,217,41DFTB 从密度泛函理论中总能的微扰展开形式出发，将总能写成带结构能量加上两体短程排斥能()totBSrepIJIJEEVR为了获得带结构能量，需构造 Hamiltonian 矩阵。onsite 对角元可以通过求解一个赝原子问题得到020()()()effIvvvrTvnr rr同时得到的赝原子电荷密度可以用来求 Hamiltonian 矩阵非对角元000effIJhTvnn，其中,IJ短程排斥势可以通过计算 DFT 自洽能量与前面得到的带结构能量的差来拟合。一种简单的方案为多项式拟合2()()0NPnnCCnrepdRRRRVRotherwise更精确的方法是使用样条函数。上面的方法实际上是 DFT 能量的零阶近似，对于电荷转移比较重要的体系不能给出较好的描述。这时，需要考虑能量的二阶修正。假设电荷密度变化可以分解到每个原子上，且具有球对称性，则修正项可写为2,12NndcoulIJIJI JEEqq。其中的对角项可以近似写18成对应原子电离势和电子亲和能之差，即 Hubbard 参数 U。U 可以通过原子能量对原子电荷的二次导数求得，或者根据 Janak 定理写成最高占据原子轨道能量对其占据数求导02atHOMOIIIIIHOMOEUqn对非对角项，如果忽略交换关联的贡献，只考虑库仑相互作用，则可写成331 88JJIIJIIJeer RrRrr，其中165IIU上述积分可以解析求得。如果采用 Guass 型原子电荷分布，同样可以解析求得通过误差函数表示的IJ。当取周期性边界条件时，须将coulE分拆为长程的1/IJR项与另外一个短程项的和，其中长程项可以通过标准的 Ewald 求和方法得到。对自旋极化可以做类似于电荷转移的处理，总能中新增一项12NIlIlIllIl I lIP P W，其中原子常数 W 同样可以通过能量对极化密度的二次导数得到，由 Janak 定理有012llIllllPWnn因为原子电荷依赖波函数，所以需要自洽求解 SCCDFTB。通常用 Mulliken 分析来求原子电荷和极化密度0occNIiiiiiiIiIqn c cn c cSqoccNIliiiiiiiIPn c cn c cS对应的 Hamiltonian 可以通过对总能变分得到，为0 1122NIKJKKIllIl lIlKlIhhsqsWWP，,IJ其中正负号分别对应自旋向上和自旋向下电子。SCCDFTB 的总能表示为nsccsccrepEEE，其中*0nsccnnnnEfc c h，不考虑自旋极化时scccoulEE。19第二章 电子结构计算方法 第二章 电子结构计算方法 2.1 原子基组 2.1 原子基组 计算机无法表示连续的物理量，所以数值计算的第一步通常是离散化。最常见的离散方法是通过实空间格点来实现，但电子结构计算用的更多的对单电子轨道进行离散化的方法是将其表示为有限个解析函数的线性组合：ia，其中被称为基函数，其集合称为基组。一种常用的方案是将基函数的中心选在分子中的各个原子上，模拟孤立原子的轨道（s,p,d,），即通过原子轨道的线性组合(LCAO)来得到分子轨道。原子轨道一般可写成(,)()(,)nlmrR r Y 其中 Ylm为球谐函数，而径向部分 Rn(r)可以选用 Slater 型函数，Gaussian 型函数，或数值函数【输入参数】。1.Slatertypeorbital(STO)STO 是类氢原子轨道，所以直观上是最佳的基函数。其表达式为 1,nrlmrreY 其中 n 为主量子数，相当于把类氢离子径向波函数正交化限制去掉，而要求其指前因子为单项式，以简化计算。Slater 型函数是非常好的基函数，但是计算三中心和四中心双电子积分比较困难。ADF 程序使用 Slater 型基组,其多中心积分都是采用数值积分方法。另外，使用Slater 型轨道只需要确定 Slater 指数就可以得到一个基函数，所以对重金属元素很容易得到可用的基组，因此在重金属计算中常常被使用。2.Gaussiantypeorbital(GTO)我们看到，Slater 型轨道的缺陷是 HartreeFock 方法中最经常出现的多中心积分无法解析求解，那么什么函数才能胜任呢？1950 年，Boys 提出用 Gaussian 型函数作为基函数来展开分子轨道。通过在 Gaussian 函数前面加上不同的因子，可以得到对应于 s、px、py、pz、dxy、dyz、dzx 等的 Gaussian 函数。2(,)ijkrx y zx y z e i，j，k 表示轨道类型，如下表所示：i j k i j k =0 0 0 0 s=2 1 1 0 dxy=1 1 0 0 px 1 0 1 dxz 0 1 0 Py 0 1 1 dyz 0 0 1 pz 2 0 0 dx2 0 2 0 dy2 0 0 2 dz2 因为 x2+y2+z2是球形的，六个 d 轨道可以表示成 5 个正则 d 轨道：xy,yz,xz,x2-y2,3z2-r2。GTO 乘积定理：双中心函数乘积等于单中心函数（两个高斯函数的积仍是高斯函数）。20对于两个中心在 A 和 B 的 Gaussian 函数，其乘积为中心在 P 处的 Gaussian 函数，P 到 A 和 B的距离由 Gaussian 函数的指数系数确定，其指前因子是由指数系数和 AB 间的距离确定。这样通过 Gaussian 函数可以把双中心积分变成单中心积分。相对于 Slater 型函数，Gaussian 函数还有一个优点，Slater 型轨道的指数上是 r 的函数，而 r 等于(x2+y2+z2)1/2，这样就把 x，y 和 z 三个方向耦合起来了。而 Gaussian 函数是指数为r2的函数，从而可以分解成 x2+y2+z2，这样就可以把三维积分变成独立的三个一维积分，可以大大简化计算。相比与 Slater 型轨道，Gaussian 函数在 r 大的时候衰减太快,而在原子核处又没有尖点(cusp)。所以，GTO 严格说来并不是一种轨道。另外，高斯函数没有节点，所以肯定没法只用一个高斯函数来表示价轨道。3.收缩高斯基组 STOnG 为了弥补 Gaussian 函数的缺陷，我们可以用几个指数不同的 Gaussian 型函数拟合一个Slater 型轨道(在电子结构计算中固定展开系数，所以称为 Contracted 高斯基)，这样在 r 大的情况下的衰减特征基本上能趋近原子轨道。若用了 n 个 GTO，则记为 STOnG。那些被用来拟合 STO 的单个 GTO 就叫原始（Primitive）基函数：,;,NGTOaaaCX Y Zi j k 一个常用的小基组是 STO3G，但是 STO3G 通常不具有足够的灵活性来描述价轨道的电子分布。4.分裂基组与分裂价基 splitvalencemulti 如果将一个原子轨道用几个 STO 来表示，则称为分裂基组。例如，双基组表示每个原子轨道用两个 STO。通常，内层轨道可以通过一个 STO 来获得较好的描述，而价轨道偏离原子轨道较多，往往需要用几个 STO 来表示。这样，我们可以只对价轨道进行分裂，称为分裂价基。如：321G 是 valencedouble基组，其中内层轨道是由 3 个 contractedGTO 来拟合的一个 STO，而外层的价层轨道则是由两个 STO 基函数来描述，它们又分别由 2 个和 1 个GTO 来拟合。5.极化基组 Polarizationfunction 对处于各向异性环境中的原子，有时需要使用极化基组。例如：如果不加极化函数，计算得到的 NH3的结构是平面型的，但实际上是三角锥结构。对 H2O，dxy轨道可以显著加大重叠：21631G(d)(也写作 631G*)对非氢原子添加 d 函数。631G(d,p)(也写作 631G*)对氢原子添加 p 函数且对非氢原子添加 d 函数。6311G(3df,3pd)在非氢原子上加的三个 d 函数和一个 f 函数，在氢原子上加的三个 p 函数和一个 d 函数。6.弥散函数 Diffusionfunction 弥散函数是在空间比较弥散的 s 和 p 轨道的基函数(把指数很小的函数添加到基组中。应用:含有孤对电子的体系，阴离子（Anion），以及其他带有明显负电荷的体系，激发态体系（excitedstate），含有低的离子化能的体系，以及纯酸的体系，计算极化等。631+G 对非氢原子添加弥散 s 和 p 函数 631+G 对氢原子添加弥散 s 函数且对非氢原子添加弥散 s 和 p 函数【作业：用 6311+G(3df,2pd)计算 H2O 时，用到多少个轨道，多少个 GTO？】7.关联一致型基函数 ccpVDZ，其中 cc 为 correlationconsistent，变分优化时不限于 HF but calculations withcorrelation。在关联基组的基础上也可以添加弥散函数,如 augccpVDZ:one set of diffusefunctionforeachangularmomentum8.数值原子基组 NumericalAtomicBasisSet 直接作（赝）原子轨道计算，存储径向函数，样条插值。优点：分子解离成原子 exactly。ConfiningPotential,energyshift,局域基组。Dmol3 软件包：采用双数值基(DN)、双数值基加极化扩展轨道的双值扩展基(DNP)、双数值基加 d 类扩展轨道的双值扩展基(DND)。Siesta 软件包：采用 DZP 基组、TZP 基组等。极化、扩展等轨道可以通过离子、激发态、类氢原子、电场微扰等方法产生。2.2 平面波基组 2.2 平面波基组 任意单电子波函数都可以写成平面波叠加的形式：()innCedg rrgg 当 g=0 时波函数为常数，g 越大时振荡越厉害。将平面波展开离散化往往需要借助于周期性边界条件(PBC)。晶体本身具有周期性边界条件，而对分子(1D)和表面(2D)体系，可以分别通过 supercell 和 slab 模型来施加周期性边界条件。对 PBC 体系，根据 Bloch 定理，离散化可以分两步走。2.2.1 单电子波函数的离散化 2.2.1 单电子波函数的离散化 Step 1：Bloch 定理定理 势场具有晶格周期性时，波函数可表示为 ()innuek rkkrr 其中)(rknu为具有晶格周期性的函数。合理性：因为 222hvm r 22所以晶格平移算符和哈密顿量可对易：,0t h 111()nnnt rRrr 因为 2121nnnn，所以 的一个可能形式为exp()nik R。根据 Bloch 定理，单电子波函数有一个周期部分，其傅立叶展开可以离散化（倒格子）。于是，整个平面波展开可写为 ,innCek G rkk GGr 其中 G 为倒格矢 b b 的整数倍。因为2ijija b，k 只需在第一布里渊区取值（周期R2）。在实际计算中需对G进行截断，而平面波基组的一个优点就是可以通过增加截断能量（cutoffenergy）【输入参数：ENCUT】来系统地改善基函数集的性质。Bloch 定理也可以看成:将波函数按 k 分类。一般说来，对每个 k k，占据轨道数等于单胞内电子数目。于是，求电荷密度由对无穷多占据轨道求和转化为对布里渊区积分:*occnnnBZd d kkrrr k Step 2：Born-von Karman 边界条件边界条件 引入 BornvonKarman（波恩冯卡曼）边界条件可以进一步离散 k。考虑一维情形 假设存在一个足够大的周期 N，这时，一维链可以看成一个环 21ikNRnN nekmNR，其中 m 为整数,取值范围为(,)22N N或(0,1)N。这样，在倒空间第一布里渊区里 k 取分立值了。N 越大，k 越密。上述两步离散法对应着周期性体系中两组对易关系。第一组（蓝色）：cutG对应着实空间格点的疏密，二者互为傅里叶变换。对动能相关项，在倒空间计算比较方便，而势能相关项通常在实空间处理比较方便。所以，在计算哈密顿量矩阵元的时候，经常需要做傅立叶变换。另一组对易关系（红色）:实空间中的 BornvonKarman 边界条件对应着倒空间中布里23渊区内 k 点的疏密，从而决定了布里渊区积分的精度。2.2.2 Wrap around error 2.2.2 Wrap around error 假设波函数截断波矢为 Gm，现在考虑其他物理量所需的最大傅立叶展开分量。因为电荷密度在倒空间为两个波函数的卷积，将包含最高达 2Gm的傅立叶分量（对超软赝势需要更高傅立叶分量）。对有效势 veff，我们无法估计交换关联势的最高傅立叶分量，但是只有最高到 2Gm的有效势分量是有意义的：,|()neffeffnnvv GG GGGG同样，由卷积定理，哈密顿量作用在波函数上时，由于有效势有 2Gm分量，波函数有 Gm分量，乘积将有 3Gm分量。取 3Gm的 FFT 网格时可完全避免 wraparound error(WAE)。现在的问题是，如果只需确保最后得到的乘积波函数的最高到 Gm分量的傅立叶系数是正确的，此时傅立叶变换的网格需要多密【输入参数：NGX/NGY/NGZ】呢？分析表明，当 FFT 网格不够大时，WAE 产生的误差有一个折叠效应（参见附录）。这样，取 2Gm的 FFT 网格时，WAE 会影响到乘积中所有大于 1G 的系数，但不影响 1G 以内的部分，所以，不影响到波函数。2.2.3 布里渊区积分布里渊区积分对周期性体系，电子结构计算与分析都涉及到布里渊区(BZ)积分，其一般形式为：1()BZnnnBZdkkk通过波恩-冯卡曼边界条件，在实际计算中积分变成求和的形式。更一般地，积分可以变成对部分特殊选取的 k 点进行带权重的相加：1BZBZwkk24这样，就涉及到 k 点的选取。在计算过程中，要用尽量少的 k 点得到尽可能精确的积分。我们考虑如下问题：怎样对有直到 N 个傅立叶分量的函数精确求积分？以下是对 2 个傅立叶分量的积分：2012()sin()sin(2)fkAAkAk 220022()3 ()()22Ifk dkAfkfk 对任意函数，如果取如上两个 k 点，可以保证其傅立叶展开以后直到 2 阶分量的积分为零。因此，为了获得较高的效率需要选取一些特殊的 k 点，原则是对一定的 k 点数量使得尽可能高的傅立叶分量的积分都为零。同时 k 点的选取也需简便易行，综合考虑后较好的方法为Monkhorst-Pack9取点法，它将三维问题转化为三个一维问题。一套 MP 网格是 BZ 内均匀分布的规则格点集，其数目由参数1