分子轨道计算.docx

分子轨道计算第一节引言分子中电子结构理论的进展几乎是与电子计算机的发展平行的.换句话说，分子中电子结构的理论 研究在很大程度上依赖于电子计算机的应用.因此，分子轨道计算就成为计算物理学中的一个重要方面。众所周知，分于是由原于组成的，原于之间的相互作用使分子的电子结构比原于的电子结构复杂得 多.在固定核或玻恩一欧本哈默的近似下，分子的定态薛定谔方程为(1.1)其中p、q标记原子核，i、k标记电子；Rpq为核p和核q间的距离rik为i电子和k电子之间的距离， 而rqi为核p和i电子之间的距离如果长度和能量都采用原子单位，方程(1.1)可以简化为(1.2)相应的哈密顿算符为。S 态 E 声 E (1.3)方程式(1.2)可写成*财=£叽(1.4)方程式(1.4)具有能谱:(1.5)各能级对应的波函数分别为W。，W，W2，其中E0为基态能量.在量子力学中常用变分法来求解薛定谔方程.这一方法的基本原理是：对任何满足一定边界条件的 尝试波函数中，成立关系式(1.6)这个关系式的证明是容易的.我们知道，哈密顿算符的本征函数构成一个正交归一的完全集合，亦即甲回 >5=1(1,)将H算符对其本征函数集合展开，我们有(1.8)上式两边同时减去E0可以得到£*0 =平(£矶)| 皿叭| 可*一矶在中态的平均值为9“7.|。=密*矶)皿咐，.(1.9)因为E0是基态能量，从(1.5)可见对一切k有Ek-E0.所以(1.9)式的右端不会小于零，即Fo|d»>o亦即这就是(1.6)式.根据变分原理，对于任何一个含有参量气，入2,，勾的尝试波函数中，(1.6)都成立.因此，我们 可以选择参量气，入2,，与，使|H|/|取极小.换句话说，要求这些参数满足条件|-?=0耳=1,2,m.(1.10)从方程组(1.10)解出气，匕，勾，就可以得到近似的波函数.用这个方法，相当于先求出最低的基 态能量E0和基态波函数 再利用尝试波函数中与W°正交:|中0；= 0的要求，可以得到激发态的E1依此类推，可得到一系列的本征值和本征函数。通常所用的尝试波函数具有形式O = G0+C,+Cn<Pnr (111)(1,11)其中和是巳知的选定的基函数，它们应满足波函数的边界条件.而这些基函数的线性组合的系数就 i，Cn.于是相应的能量是尝试波函数含有的参量气，妇，勾.通常将这些系数记为C，C2,IE- 9/IS 一 部:GH"w(1.12)其中短阵元Hti= «2><|jr| 0/> ；S</= <0><| 0/>,将能量E对系数C求极小，得到一组n个久期方程 1£ (Zf 一ES)Cj=0(:=1,2, »«),(1.13)上式仅当系数行列式为0时才有非零解，即det(HfL£Sv) = 0. du)(1.14)它的本征值E和本征向量C 就构成所求的薛定语方程的能量和波函数第二节Roothaan方程和从头计算对于多电子体系，直接从(1.4)出发利用变分法求解，实际上是很困难的，必须再作进一步近似.在量子 力学中，利用单电子近似；可从(1.4)得到HartreeFock方程尸物=&+ 客们=田礼(« 1,2, ,«),(21)其中库仑算符和交换算符分别为 以1)=/(2)印(2) - dr3；(2.3)Kj(l 油(1)=J 牝(2)奶(2)=仞 d。此处及以下dT2均为三维空间积分元，为简单起见只用一重积分符号表示；函数中的宗量1, 2分别表示和r2.而H算符是(2.4)算符(2.2)和(2.3)中都含有未知函数.实际上，(2.1)是三维空间中的非线性积分微分方程.一般都要采取 选代法，也就是物理上常称为的自治场方法来求解.现在对(2.1)采用前节的变分法来求解.这首先要求选择了组恰当的基函数.因为分子是由原子组成 的，所以很自然地采用原子轨道为基，将分子轨道用一组原子轨道的线性组合来表示.设原子轨道为，说，队，构成一组基集合.将分子轨道表示为京勺线性组合(2.5)当然，系数C .应满足使中归一化的条件.根据上节，从变分法得到久期方程 卬iTrt，1(F“一&SQ = O,(2.6)这里£ i为单电子的轨道能量，也称为分子的轨道能级，矩阵元(2.7)方程(2，6)式是齐次的线性方程组，称为Roothaan方程，它有非零解的条件是(2.8)这是标准的短阵的广义本征值的问题，它的求解有标准于程序可以调用.困难在于矩阵元尸"的计算。则在闭壳层的情况下，分子含有2n个电子，经过繁复的为简单计，假设原子轨道都是实函数组，解析推导，我们得到旦尸:壹Pg (如如)土 (以"),(2.9)其中与电子之间相互作用有关的量(如1人。)="伊“(1)啊衣(2)代(2)打0弓(2.10)为与原子轨道有关的六重积分，而!(2.11)是所谓密度矩阵.它是从电荷密度表达式和(2.5)得来。由电荷密度表达式可以得到p(r)=2 习耐(r)=2)P“ 构(丁)。”(>*),(2.12)因为pG )在全空间的积分等于分子的总电子数，即2n=J“ J”得到新的P(2.13)(2.13)即构成Mulliken所谓的集居数分析，也就是密度短阵的来源。(2.9)式中的H "为IH”=侦"1)反佐d”(2.14)而算符由(2.4 )式表示.可见，在选定基函数组以后可以计算(2.10)型的积分，从而求出F"，但P四中含有未知的系数C 必须在给定一组初值的条件下才能算出F" .从而利用矩阵对角化的办法求出£ ,和C皿 和V，再进行下一轮计算.在开壳层的情况下，要考虑自旋波函数.玲和P是两种不同的自旋波函数.若体系有p个a电子， q个P电子，且p>q.于是可构成两种对应于以自旋和P自旋的分子轨道如下：#=甲矣,(2.15)(2.16)记相应的电荷密度为P，(r)=幸样(r)'=2J%3(r)%(r),<uv棉(r)V=£P么 知(r)们(r).(217)自旋密度矩阵为珞=加6,性=力C%(218)完整的密度矩阵为(2.19)于是，利用变分法原理可以分别推导出两组联立方程*(F 为一)C®<=02XF：玲 S”)C£s= 0(2.20)其中'万*=善P人。)一Pjr (四口 I 和)队,=日“+习&31和)一环(如和)./、(2.21)求解这两个方程组同样要用自治场的方法，先给定任和依的一个初始估计，计算伫和写*，然后 分别求解(2.20冲的两组方程，求出罕和£ P及波函数C*和C，通过(2.18)、(2.19)和(2.21)可以准备 新的一轮迭代.与此同时，还可以计算分子的总能量，也就是在单电子近似下，§2.1中的(1.12)所表示的能量.在 闭壳层情况下，分子的电子总能量为IE=S)P”H“+g 2)巳/(“胡|如) Z “AO*(2.22而在开壳层时，分子总能量为I(2.23)E=gF"H“+；翠P“PkSI和)-(P：0 +PFL)(gMy)这种直接计算、和，"，然后求解Roothaan方程的方法，叫做从头计算法(26 initio方法).这个 方法的主要优点是，它把原来求解一个积分微分方程(2.1)的问题化为一个代数方程(2.6)的求解，使分子 结构的计算迈进了一大步.但是这个方法的计算工作量非常大.计算量主要在于对(2.10)型积分的计 算.这是一个六重积分，它的被积函数可能牵涉到四个不同的原于轨道Pp，甲v，平人和中°，而这四个 轨道可能分属于四个不同的原子，这就构成所谓多中心的积分.如果基函数组包含m个原子轨道，一 般讲来，可能有多达m4个(2.10)型的多中心多重积分需要计算，其工作量可想而知.因此，这个问题成 为分子结构计算中著名的难题.尽管在50年代初期已经得到Roothaan方程(2.6)，但实际应用它来求解 分子的电子结构在当时简直是望尘莫及的.只有出现了相当大型的电子计算机以后，从头计算方法才有 可能实现.在具体的求解过程中还遇到原于轨道的选择问题，一般采用原子轨道的计算.用那些方法算出的原 于轨道可以作为基函数.这样的基函数都是离散的数值表.因此，(2.10)积分和交迭积分(2.7)的计算都 必须用数值积分形式来完成.这样的方式有时又叫做离散的变分法，它完全是建立在离散化的基础上的.另一种考虑是用解析形式的基函数来模拟原子轨道.这就可以避免数值积分.对(2.10)作一系列的 解析推导，可以适当减少计算工作量.常用的解析形式的原子轨道有Slater轨道和高斯(Gauss)轨道两种.1. Slater 轨道Slater轨道的选取形式为如-即出(静)(2.24)的函数，这里A是归一化常数*=(2£)平 C(2”)!T, n、l、m依次为主量子数、角量子数和磁量子数，Z称轨道指数.可以用拟合的办法选取Z，使(2.24)更 接近真正的原子轨道，下面是对简单的双原子分子氮化氢HF的计算实例，计算时将F原子核固定在原点，H位于(0，0，R)，R是键长1.733au或0.9171A这里au表示原子单位，即玻尔半径a。.选取基集合为F原子的1s，2s，2px，2py和2pz轨道，以及H原子的1s轨道.这 是最小基集合轨道.采用的Slater函数为'''对F：对H：而Z = 8.7； Z2=2.62；，3 = 1.0。这是闭壳层分子，总共10个电子，共占有5个轨道.6个基函数得到6 阶的久期方程组.经过6次迭代收敛.在迭代过程中总能量的变化如下：第一次迭代：-103.93973au第二次迭代：-104.66070au第三次选代：-104.67676au第四次迭代：-104.67184au第五次迭代：-104.67184au第六次迭代：-104.67184au最终得到的分子轨道总能量和分子波函数系数C见表2.2-1W衰 2.2-1q (at )原子轨卜、-26.139-1.4760.566-0.465'0.46S-0.47710.9963-0,24350.0839000.080020.01630.93220.471500二0.669930.00240.09070.687000-0.806540001,0000、o60001.000 L06-0.00460.16060：576101.05022.高斯轨道高斯轨道与S1later轨道的差别在于其径向部分的指数衰减呈r2的形式，即为 pn-1e -a r 2的形式.角向部分仍为球谐函数Yjp,甲).可见，随着r的增加高斯轨道的衰减比Slater轨道快得多， 这是有利于多中心积分的计算的.而且选用高斯轨道可以将形如(2.10)和(2.14)的多中心积分全部用解析 形式推导出来，可使计算大为简化。因此，在多原子分子轨道的计算中，人们往往来用高斯轨道为基， 以减少计算工作量.从球谐函数的表达式中可以看出，在直角坐标系中 爪1丫而使甲)实际上是一个多项式的形式，因此 往往将以原于A为中心的Gauss轨道函数定义为尤(4 Ct, r, m 矽=Nx 我 2% e ”。(2 25)2t(l+«+n)+»/a a：+m+n+l/l其中N是归一化因子；Ir2I(I+w+n)+s/s a+w+n+i/i-ii/s(2.26)实际上，Slater型的函数(2.24)式能较好地模拟原子轨道，但用Slater函数作基的多中心积分的计算非 常复杂.为弥补高斯轨道不像”原子“轨道的缺点，往往增加基函数组的数目.例如，先用几个高斯轨 道来拟合32ater钦道，也就是将Slater函数用高斯函数展开的办法.利用高斯轨道可以将(2.9冲的多中心积分(2.10)和(2.14)全部用解析函数的形式表示出来，但其计算 仍是非常复杂的，只能用电子计算机来完成.所以从头计算法在很大程度上是依赖于电子计算机的.为了减少计算工作量，人们考虑了很多节约的办法.一个常用的办法是所谓冻结芯态.因为当原子 组合成分子时，其芯态电子是变化很少的，活跃的是价态电子，所以常常将芯态电子“冻结”起来，只 考虑价态电子的运动，这样可以减少基函数组的数目.第三节几种常用的近似方法因为从头计算法的工作量太大，特另U对大分子的计算非常困难，所以六十年代以来发展了分子 结构计算的近似方法.下面介绍常用的几种近似方法。一、CND0方法CND0方法的全称是完全忽略微分重叠(Complete Neglect of Different Overlap)方法.这种方法首先只处理原子的价电子，将芯态电子完全冻结，基函数组就取为价电子的Slater轨道.这种近似的基本思 想是将不同原子轨道乘积的微分看成零，即9?Fdr=0（“与助）.于是在久期方程（2.6）中的置叠积分就不出现了，因为S四=8即。方程（2. 8）化为det（F 好一眼”）=0,（3.1） 而且电子排斥积分简化为（如如）=（加办）8“北，（3.2） 上式表示四中心积分披简化为两中心积分进一步再假设非零双电子积分只与*和也所属的原于A 和 B有关这样就可以用A 和 B的s轨道来近似*和也，使计算过程更为简单此时（3.3）:呻= y ab是依赖于原子A和B的s轨道的积分.至于（2.9）中的H"可引进下面的近似，设p轨道属于原子A，则HZfSf （3.4）*）（3.5）是仅与原子A有关的量，其中没有相互作用量.记乙，=（，1乙也），（3.6） 根据上面的原则，这里的%轨道可用A原子上的s轨道计算之若M和v是同一原子A上的不同轨道，则根据忽略微分置叠的原则，H叫=0（当日。V，而同属原于A时）（3.7）如果P和v分属不同的原子A和吓，这时不能再忽略,而假设它眼原于和B原于之问的重叠积分、 成正比，即乩尸服疝，（3.8）P人6是一个适当的参数，它标志原子A和B之间的成键关系，称为成键参数，综上所述,可将期表示的矩阵元七简化为（3.9） § p 冲）n,+ 浴'死斤妣 一v 尸“=们舟厂§“皿 （件*V）,而U州是原子参量,咏是成键参量:(3.10)上式表示求和取遍原子龙A上的所有轨道.相应地，总能量的计算公式为(3.11)以上各式也可以类似地推广到开壳层的情况，这里不再赘述了。从上面的计算公式可以看出，因为P"中含有未知系数Cp，所以整个计算过程仍需要用自治的 方法来完成.其次公式中含有参量U辿和P A3,这些量往往是用经验的方法来选取。例如，将U辿取为原子的电离势，或者取为电离势I p与亲合势Ap的平均”叩=顷+皿(312).更有将VAB取血=2砰m (3.13)根据对这些经验参量的不向选取的方法，又使CND0方法区分为CND0 / 1，CNDO/2 等等.下面是用CND0 / 2方法计算NH2分子的平衡键角的实例.基函数组取为H的ls轨道，N的2s，2p轨道，共计6个函数.巳知N、H之间的键长为1.024A.分 子NH2共有七个电子，是开壳层分子，可能有两种组态，NH2属于Cm点群的对称性，故分子轨道的对 称性为和ab1和a2、b2.分子可能的两种组态，一是在a1上少一个电子，二是在b1上少一个电子， 将前者记为A1，后者记为B1.计算得到这两种态的总能量随键角。的变化关系如图3.1所示.-awtrK)圆堤图3.1 NH2两种组态的总能量随键角。的变化 当然.还可以求总能量与键长的变化关系，那就需要计算位能面。第四节多重散射Xa方法多重散射Xa自治场方法(MSXaSCF)又称散射波Xa方法(SWXaSCF).这种方法和前面选 取原子轨道为基，以变分法为基础的方法是不同的。它采用了固体物理能带计算中的很多想法，形成与 原来量子化学不同的新方法.这种方法最初的想法是S1ater提出的，后来由他的学生K.H.Johnson经过多年的努力，逐步完善，并在电子计算机上具体实现，获得了成 功.下面叙述xo方法的具体内容.首先，Slater对HartreeFock方程作了修改，将交换势用自由电子气的统计平均来作近似.用Rydberg单位后(4.1)在这种修正下得到的方程就称为Hartree一FockSlater方程= E 虬(4 2)(4.2)+(4 3)(4.3)(4.4)V0(r )是库仑势(4.5)颌)是电荷密度:服)=j碧阵Vex(r )由(4.1)式表示.(4.6)ni是第i个分子轨道上的电子占有数.1 Xo方法采取重要的简化措施是对势场V(r)作所谓Muffin-tin近似.将全空间划分为三个部分(见图 4.1)区域I是由组成分子的各个原子的近邻(称为原子球)构成，这些原子球相互接触(相切).然后用一个 大球(称为Watson球)将这些小球包围起来.在大球外部就是区域II，余下的球间部分(即大球内部减去 区域I的部分)为区域II可见区域I和II都是球形区域.所谓Muffin-tin近似就是在这些球形区域中对势场作球面平均近似，而 在区域II中作体平均近似.设区域I中的原于球用j标记，区域II用0标记，在这些球中势场的球面平 均为匕(门)=£(及+r,)dS(4.7)R.是第j个原子核的位置，是以第j个原子球心为中心的局部坐标.在区域II中就取势场为一常数U:%)(4.8 )在以上的势场Muffin-tin近似下，我们可以将方程(4.2)的解分区用多重散射波的形式表达出来.首先，在I区和II区，由于势场仅与径向分量r.有关，故可将方程(4.2)式分离变量，得到形式解(4.9)PS)CM(E)匕 3, 其中R(E,)是径向方程-寿寿何器）+ - +服力-归传。(4.10)的解，（4.9）和（4.10）中的j = 0就表示II区的解，即儡（儿）=:SC焉用（尻r。泪（*）.（4.11）这里,表示r的角部分.在区域II中，由于势场是常数，所以分离变量后的径向方程可以用贝塞尔函数 解出来.于是，记L=（l，m），贝0争再思金（奶）儿（£）+*避心侬r°）YL（£）当E>刀税（研）印（如。）"（上）当 H,J 匕上（4.12）其中k 二幅-E1/2，ji（x）是球贝塞尔函数，ngx）是球诺伊曼函数，即*）=/乎妇 G），m）=J专叫，砰）（”）=/专右边是半整数贝塞尔方程的解，分别为贝塞尔函数，诺伊曼函数和第一类汉克尔函数.而IAG）=iT，Kbf）,研（）= "矽）（技）,区域II中的解（4.12）具有散射波的形式，它是各个原于球的散射波和大球（III区）散射波的送加.但另一 方面，如果从第j个原子球来看，区域II的解又是第j个球的入射波和散射波的和所以，在区域II中的解又可表示为（当j。0时）V»n（r） =老5功（奶）L伤）+孕4挪）（研）"（A）当E<yu顷"g（泌呼杯河（£）当£>九（4.13）上式中第一项是入射波，第二项是散射波.（4.12 ）和（4.13）所表示的函数应该是一致的，从而经过解析推导，可以得到关系式I班'=定 G 总（E）成;+gS£”（4.14）当j = 0时（即对外球而言），相应于（4.13）的【1区波函数|3跳虾（她比航。）+、刀出（研）皿（；。）当E<yn饥i（r）= L入 '八孚更M研）皿。）+羊盗讯刑。比乙（上）当E泓,故相应地有I班船（E）理,.（4.15）当E 匕时以寮E)=U 8"旧 T)je咨玲仁”嫩 ) 睥如),S 纷(旧)=4七)(一 1)+工人”(£,)标以R，o)VzX£。)，(4.17)S%(E) = S%(E) s)而当E 匕时I朝£)=勒)4浦编哪所，)膈(疆,(4.19)IS技(E)=4"-" S日叮切(L, )力，(利殆。)丫匕扁),(4.20)sy=s“E), (4.21)其中，L)= |k(r)耻，(r)*/(r)dr(4.22)这样得到的三个区域中的分块解在区域的边界上应该相互衔接，这是置子力学的基本要求，即波函数应 该是连续可导的。我们要求不同区域的波函数在跨过区域边界时函数本身及其一阶导函数应该连续.根 据这个要求，可以得到系数之间的关系式C 见(E,S)I 挪)(物)+ 助)(423)(4.23)即AABME),恐邮)=斑.(4.24)止_ 3(知)，死(E0).9 一脸饱)，死(E0)r(4 25)(4.25)以_ 研g顼?(E,6。)讯也),&?(£&) ' (4.26)当E 匕时你£)=一场(砌)，昭(E,b况3(他)，州互方)(427)(f匕/)40( P”(泌0),殆(£,如)I3(他)顼谯,儒)'(4.28)这里b.为第j个原子球的半径，b0大球半径，还用了朗斯基符号【甲)=冲一/虬(4 g)(4.29)将(4.24)代入(4.14)和(4.15),得到关于系数Aj的线性关系登夺(£)纱思，一善S%(E)/ = o；岑弟s”(E)左，一善口?3)冲心理，=o.它应该有非零解,于是系数行列式应为0,即L(E)纷一 S/,(E)S 螃(E)挪(E)T&“(4.30)=0,(4.31)其中丑务=<3功宓.(4.32)方程(4.31)就是Xa方法的久期方程.它形式上很简单，但和普通的矩阵本征值求解有本质的差别，矩 阵的每个矩阵元中都含有E,并且不是E的线性函数，而是隐含有E的特殊函数.因此普通关于矩阵 对角化的办法在这里不适用，求解极其复杂，这是Xa方法的一个难处.这些未知的能量互都是通过径 向方程(4.10)的解表示出来的.一般在用Xa方法求解时，都要先利用分子的对称性，对久期方程(4.31)作对称性约化.这样不仅 可以降低矩阵的阶数，而且避免了重根，求根也就容易一些.但又因势函数中台有未知的波函数，所以 求解以求解需要通过自洽的过程.Xa方法的程序比从头计算方法复杂，但计算启则减少很多。第五节固体电子结构计算固体电子结构的计算比分子的情况更要复杂得多.对完整的晶体，由于原于排列的周期性，能带理 论的计算现在已经发展得比较完善.如利用原子轨道的Bloch波为基的紧束缚方法(LCAO方法)，当原 子轨道取为高斯型函数时往往又称为LCGO方法.从前一节关于分子结构的计算中，可以看出这种方 法的计算工作量很大.于是又出现了从头计算的和经验的紧束缚两种方法.后者往往根据已知的实验值 采拟合一些参数.在一些复杂问题的计算中还是很奏效的.在固体电子结构计算中，还有采用以平面波为基的方法，也就是利用傅里叶展开的办法.在此基础 上发展起来的有OPW(正交化平面波)方法和赝势法，以及后来取得很大成功的自赝势方法。此外，还有在对势场作Muffin-tin近似下导出的APW(增强平面波)方法和利用格林函数的 KKR(Korringa，Koho，Rostoker)方法.这些方法的主要困难在于久期方程是非线性的，类似前节由Xa 方法导出的久期方法，求解跟困难.后来提出了线性化久期矩阵的想法.这就是所谓LAPW(线性化增 强平面波)方法，以及FLAPW(全势场线性化增强平面波)方法。后者可以说已经丢弃了 Muffintin近 似.当然，这些方法的计算工作量很巨大，没有大型计算机是不可能进行的.关于这些方法都有专门的 文献可查，这里不再叙述了.关于完整晶体的电子结构已有较好的研究，近年来比较感兴越的则是不完整晶体的电子结构.例如， 表面和界面问题、杂质问题、以及非晶的电子结构问题，对于这些物体，通常的原子排列的周期性已经 被破坏了，至于固体的表面问题，尽管在平行于解理面的方向仍保持二维的周期性，但垂直方向的周期 性巳经破坏.为研究表面(包括清洁表面，以及有吸附物的表面)的电子结构，人们常采用所谓薄片(Slab) 模型(图5.1)，它是由具有二维周期性的若干原于层作为基本单元，加上若干真空层，构成一个周期.然 后沿垂直方向作人为的周期性延拓，从而得到了三维的周期性.于是能带计算的方法又可以应用于这种 薄片模型.当然现在处理的是比完整晶体大得多的原脑.计算的复杂性和困难当然也要随之增加很多。 关于薄片模型的另一种处理原子展.真空是图5.1薄片模型另一种值得一提的方法是所谓集团模型的方法.不论是表面问题或者杂质问题所引起的电子态的变 化都是相当局域的，而人们最感兴趣的也正是这些局域的电子态.因此，近年来有人采用有限个原子构 成的集团来模拟无限的晶体.这些原于是晶体中的一部分，它们的位置就职为晶体中的几何构型.这种 原子集团已经没有周期性，对它的电子结构，本章前几节叙述的分子轨道方法都可以应用.计算的许多 结果说明，在一定程度上这种集团模型的方法对于研究清洁表面和吸附表面的几何构型和电子结构都是 相当有效的。当然，在具体处理时，还有许多细节问题，如集团的选取方法，集团边界的处理等等， 还需要对具体问题作具体的分析研究.