第一性原理计算原理和方法.doc
- 文档编号:2504475
- 上传时间:2023-05-03
- 格式:DOC
- 页数:44
- 大小:918.50KB
第一性原理计算原理和方法.doc
《第一性原理计算原理和方法.doc》由会员分享,可在线阅读,更多相关《第一性原理计算原理和方法.doc(44页珍藏版)》请在冰点文库上搜索。
第二章计算方法及其基本原理介绍
化学反应的本质是旧键的断裂和新建的形成,参与成键原子的电子壳层重新组合是导致生成稳定多原子化学键的明显特征。
因此阐述化学键的理论应当描写电子壳层的相互作用与重排,借助求解满足适当的Schrodinger方程的波函数描写分子中电子分布的量子力学,为解决这一问题提供了一般的方法,然而,对于一些实际的体系,不引入一些近似,就不可能求解其Schrodinger方程。
这些近似使一般量子力学方程简化为现代电子计算机可以求解的方程。
这些近似和关于分子波函数的方程形成计算量子化学的数学基础。
图2-1分子体系的坐标
2.1SCF-MO方法的基本原理
分子轨道的自洽场计算方法(SCF-MO)是各种计算方法的理论基础和核心部分,因此在介绍本文计算工作所用方法之前,有必要对其关键的部分作一简要阐述。
2.1.1Schrodinger方程及一些基本近似
为了后面介绍各种具体在自洽场分子轨道(SCFMO)方法方便,这里将主要阐明用于本文量子化学计算的一些重要的基本近似,给出SCFMO方法的一些基本方程,并对这些方程作简略说明,因为在大量的文献和教材中对这些方程已有系统的推导和阐述[1-5]。
确定任何一个分子的可能稳定状态的电子结构和性质,在非相对论近似下,须求解定态Schrodinger方程
(2.1)
其中分子波函数依赖于电子和原子核的坐标,Hamilton算符包含了电子p的动能和电子p与q的静电排斥算符,
(2.2)
以及原子核的动能
(2.3)
和电子与核的相互作用及核排斥能
(2.4)
式中ZA和MA是原子核A的电荷和质量,rpq=|rp-rq|,rpA=|rp-RA|和RAB=|RA-RB|分别是电子p和q、核A和电子p及核A和B间的距离(均以原子单位表示之)。
上述分子坐标系如图2.1所示。
可以用V(R,r)代表(2.2)-(2.4)式中所有位能项之和
(2.5)
l原子单位
上述的Schrodinger方程和Hamilton算符是以原子单位表示的,这样表示的优点在于简化书写型式和避免不必要的常数重复计算。
在原子单位的表示中,长度的原子单位是Bohr半径
能量是以Hartree为单位,它定义为相距1Bohr的两个电子间的库仑排斥作用能
质量则以电子制单位表示之,即定义me=1。
lBorn-Oppenheimer近似
可以把分子的Schrodinger方程(2.1)改写为如下形式
(2.6)
由于组成分子的原子核质量比电子质量大103倍~105倍,因而分子中电子运动的速度比原子核快得多,核运动平均速度比电子小千倍,从而在求解电子运动问题时允许把电子运动独立于核运动,即认为原子核的运动不影响电子状态。
这就是求解(2.1)式的第一个近似,被称作Born-Oppenheimer近似或绝热近似。
假定分子的波函数Ψ′可以确定为电子运动和核运动波函数的乘积
(2.7)
其中Ф(R)只与核坐标有关,代入方程(2.2)有
对于通常的分子,依据Born-Oppenheimer原理有:
ÑAΨ和ÑA2Ψ都很小,同时MA≈103~105,从而上述方程中的第二项和第三项可以略去,于是
易知
也即该方程可以分离变量而成为两个方程
(2.8)
(2.9)
方程(2.8)为在某种固定核位置时电子体系运动方程,而方程(2.8)时核的运动方程。
E(R)固定核时体系的电子能量,但在核运动方程中它又是核运动的位能。
此时分子总能量用ET代表。
因此,在Born-Oppenheimer近似下,分子体系波函数为两个波函数的乘积(2.7)式。
分子中电子运动波函数Ф(R)分别由(2.8)和(2.9)式确定。
电子能量E(R)为分子的核坐标的函数,从(2.9)式看出它又是核运动的位能。
在空间画出E(R)随R的变化关系称为位能图。
l单电子近似
体系的电子与核运动分离后,计算分子的电子波函数Ψ归结为求解下面的方程
(2.10)
(2.10)式是量子化学的基本方程,目前已有多种求解这个方程的方法。
这些方法的区别首先是构成Ψ的方式及其相应的近似。
最常用的是Hartree建议的单电子近似[6]。
在多电子体系中,所有电子势相互作用的,其中任意电子运动依赖于其它电子的运动。
Hartree建议把所有电子对于每个个别电子运动的影响代换成某中有效场的作用。
于是每个电子在核电荷及其余电子有效场产生的势场中运动仅依赖于电子坐标。
从而,电子运动分开了,对于多电子体系中每个电子可以引入单电子波函数,这种单电子波函数是(2.10)式单电子Schrodinger方程的解,其中含有算符1/rpq项,用只依赖于所研究电子坐标的有效场代替。
整个多电子体系波函数等于所有电子的单电子波函数(轨道)乘积。
电子还具有自旋角动量s,其分量sx,sy和sz满足普通角动量算符的对易关系。
算符s2和sz完全给定了电子的自旋,电子自旋波函数h(x)满足方程
(2.11)
其中x是自旋坐标,通常把对应于自旋1/2的波函数记为a(x),而把自旋ms=-1/2波函数记作b(x)。
在非相对论近似下和不存在外磁场时,电子的自旋和空间坐标无关,因此,因此电子的自旋轨道可取成
(2.12)
考虑到自旋变量的多电子波函数由自旋轨道组成,他应当是体系总自旋S2及其Sz的本征函数
(2.13a)
(2.13b)
构成体系多电子波函数Ψ时,必须考虑Ψ相对于任一对电子交换的反对称性要求,此所谓Pauli原理[7]。
因此,一般不求出Hartree方法的简单乘积型波函数Ψ,而是求出对应于按自旋轨道电子的所有可能置换方式的Slater行列式波函数,此为Hartree-Fock方法。
对于置于n=N/2轨道的Ψ上的N电子体系,单电子近似下波函数Ψ写为
(2.14)
该式的Slater行列式是保证反对称性要求的唯一这类函数。
引入单电子近似便确定了波函数Ψ的形式,用它可以求解方程(2.10)。
显然在一般的情况下,Ψ应当包含(2.14)型行列式的线性组合,同时满足(2.13)式的限制。
若(2.12)式中自旋部分是单电子自旋投影算符Sz的本征值,则(2.13b)式就满足。
当分子的n个轨道每个均为自旋反平行电子对占据时(闭电子壳层),一个行列式波函数(2.14)就已满足(2.13a)和(2.13b)。
对于含有未配对的电子体系,这是做不到的,此时体系波函数是对应于各种轨道填充方式(不同组态)的Slater行列式Ψl的的线性组合
(2.15)
当适当选择行列式前系数al时,条件(2.13a)和波函数的反对称性要求均可以满足。
由于存在着电子运动的相关,不明显处理(2.10)式中1/rpq项的单电子近似,完全忽略了这种相关效应,所以,Hartree-Fock单电子近似使波函数的计算产生了误差。
l变分原理
上述单电子近似只是给出了所求解体系多电子波函数的一种形式,变分法提供了求解方程(2.10)的一种方法。
Schrodinger方程(2.10)的解对应于稳定态能量。
因此若波函数Ψ是(2.10).的解,那么对于任意微小变化dΨ,取能量平均值
(2.16)
的变分应等于零,即
(2.17)
(2.16)式中积分是对Ψ的所有变量进行的,并且已假定Ψ是归一化的,即
(2.18)
由于我们寻找对应于体系基态的波函数,总能量应当是极小值。
因此,对单电子轨道施行变分就给出这种型式波函数,能量是极小值并满足(2.17)式。
从而求得的波函数Ψ就是多电子体系基态Schrodinger方程所欲求的解。
显然,为了施行变分,波函数Ψ的型式应当充分好。
两种途径可以保证这一点:
①取展开式(2.15)是从充分多项,且固定轨道Ψ只对系数al变分;②局限于尽可能少的行列式Ψl,若有可能做到就取一个,但此时把每个ψ表成可能的简单形式。
鉴于这种选择,区分出两类广泛应用的量子化学方法,价键(VB)法和分子轨道法(MO).
在价键法中,用孤立原子的原子轨道(AO)作为单电子波函数ψ去构成Slater行列式Ψl。
原子轨道的不同选择对应于不同的行列式Ψl。
对于(2.15)施行的变分,可得到确定系数al的方程。
为了充分靠近体系的能量,必须在(2.15)式中选用足够的多项,即用多行列式波函数进行运算。
用原子轨道线性组合分子轨道(LCAOMO)法提供了另外一种选择相应于体系能量极小的多电子波函数方法。
此时,对应于分子中单电子态的分子轨道ψi写成原子轨道φμ(基函数AO)的线性组合
(2.19)
实际上,这种展开有完全合理的基础。
因为靠近某个原子的电子所受的作用基本上是由该原子产生的场引起的,所以该区域中电子波函数应当近于原子轨道。
展开该式对求解变分问题的优点是明显的。
如果(2.15)式中选用极大数目的项,那么VB法和MO法就都给出同样的能量E和波函数ψ,当然表达式不完全相同。
这种唯一性的原因很简单,因为使用LCAOMO的的每个行列式均可以展开为AO组成的一些行列式。
在一般情况下,每个MO组成的行列式应展开成AO组成的所有行列式。
因为波函数ψ应通过AO组成的行列式完全集合表达,从而,当使用完全集合时,MO法与VB法所描述的ψ就等价。
当然,不用完全及表达时,两种方法的等价性就破坏了。
在极端性况下,某种方法中可以取一个行列式,此时可以直接看到MO法的优越性。
对于MO法,允许采用单行列式表达ψ(至少对于闭壳层体系),进而,通常由一些正交分子轨道组成行列式
(2.20)
其中δij是Kronecker符号。
从而是计算大为简化,并能比VB法更简单地确定(2.19)式的方程系数。
同时,MO法的基本方程能很好的适应现代电子计算机的能力。
由于这个原因,现代的MO方法已经成为最常用的计算多电子分子的电子结构的基本方法。
2.1.2闭壳层体系的Hartree-Fock-Roothaan方程
在分子轨道范围内,对闭壳层体系,在单电子近似下,用两个自旋反平行电子填充每个分子轨道ψ,可以构成一个Slater行列式(2.14)型波函数,选择轨道(2.12)的自旋部分满足(2.11)式,则保证了(1.13b)条件。
根据变分原理,若轨道ψ使得分子能量(2.16)取极小值,就求出了所研究多电子体系方程(2.10)的解。
将波函数(2.14)代入(2.16)式,并进行一些推导[见引文1-4],可得闭壳层分子的电子能量表达式
(2.21)
此处Hii是对应于分子轨道ψi的核实Hamilton量Hcore
(1)的单电子矩阵元
(2.22)
而Hcore包含电子动能算符和分子中原子核对电子的吸引能算符
(2.23)
下面两式分别表示库仑积分Jij和交换积分Kij
(2.24)
(2.25)
积分取遍电子1和2的全部空间坐标。
从(2.22)至(2.25)可以看出(2.21)式中各项的物理意义。
显然,单电子积分Hii表示在核势场中分子轨道上电子能量,由于每个ψi轨道上占据两个电子,所以乘以2。
双电子库仑积分Jij表示ψi和ψj轨道上两个电子间平均排斥作用。
由于波函数的反对称性要求出现了交换积分Kij(在Hartree方法中不考虑它),减小了不同轨道ψi和ψj上平行自旋电子间相互作用,正是它们描写了相同自旋电子运动的交换相关。
然而,在Hartree-Fock方法中,还是没有考虑反平行自旋电子间库仑排斥引起的电子相关效应。
为了求解Ψ的最优近似,必须选择一定形式的分子轨道ψi是总能量最小。
这些分子轨道相互正交,在LCAO近似下表成原子轨道ψμ的展开(2.19)式。
Roothaan最先解决了这一问题[8]。
关于ψi对AO基展开系数的方程成为Hartree-Fock-Roothaan方程(间记为HFR方程)。
下面简要推导这个方程。
首先从写ψi的LCAO展开式
分子轨道正交归一化条件给出对MO系数的附加限制
(2.26)
其中Sμν是原子轨道ψμ和ψν间的重叠积分
(2.27)
我们不难用AO基写出(2.23)至(2.25)式
(2.28)
(2.29)
(2.30)
其中Hμν是核实Hamilton量(2.23)相对原子轨道ψμ和ψν的矩阵元。
双电子相互作用积分<μν|λσ>表达式为
(2.31)
它表示电子云分布ψμψν和ψλψσ间相互作用。
从而能量E的表达式(2.21)变为
(2.32)
引入原子轨道基的电子密度矩阵元
(2.33)
上式化为
(2.34)
或者令
(2.35)
则可以用用矩阵形式写为
(2.36)
其中电子相互作用矩阵G定义为
(2.37)
矩阵J(R)描写库仑相互作用,而K(R)描写电子交换相互作用,其矩阵元分别为
(2.38)
(2.39)
作用于两个矩阵的运算Sp是指把此两个矩阵所响应矩阵元的乘积求和。
在LCAO分子轨道法中,通常假定基原子轨道是固定的,而型式不变,因此对分子轨道变分归结为对展开系数cμν的变分
(2.40)
应当指出,(2.40)式中的ψμ不变性在原则上是不必要的,已有把ψμ看成可变的一类计算方法,然而是极其困难和复杂的,常用的MO法一般不做这类计算。
当对分子轨道变分时,借助Lagrange乘法可以是能量极小化;这是极小化泛函
即变化系数cμi,使稳定点δG=0。
其中E由方程(2.34)定义。
取G的变分,则有
(2.41)
由于δcμi*的任意性,所以必须有
(2.42)
因为分子轨道在酉变换下的确定性,可以自由选择非对角Lagrange乘数为零。
定义Fock矩阵元
(2.43a)
或记为矩阵形式
(2.43b)
其中H为Hartree-Fock矩阵的单电子部分,而G为双电子部分。
于是方程(2.42)可以写为
(2.44a)
或记为矩阵形式
(2.44b)
此为HFR方程,其中E是Hartree-Fock算符本征值组成的对角矩阵,S是基原子轨道的重叠积分矩阵。
解(2.44a)或(2.44b),就给出按原子轨道展开的分子轨道系数和单电子分子轨道能量εi。
然而,上述HFR方程是数学上广义特征值问题。
为了求解方便,先借助于下述变换将它化成标准特征值问题。
由于S是Hermite矩阵,可以用酉变换θ使其对角化,即
(2.45)
从而可由其对角元平方根倒数构成矩阵S-1/2,同时S1/2S-1/2和S-1/2SS-1/2=1成立。
此时可做变换
(2.46)
(2.47)
则方程(2.44)化成
亦即
(2.48)
这是标准特征值问题,可以采用数学上标准的对角化方法处理。
显然这些方程都是系数cμi三次联立方程组,必须用迭代法求解:
当F(0)=H时,解上述方程(2.44a)或(2.44b)就得到MO系数的零级近似C(0)。
用C(0)计算储F
(1);再代入(2.44)方程确定新的系数C
(1),然后再计算F
(2)等等。
最后,当前后两次迭代所得系数C(n)与C(n-1)符合收敛精度时,迭代过程收敛。
应用现代快速计算机易于实现这样的计算过程。
2.1.3开壳层体系的非限制性Hartree-Fock方法
对于含有奇数电子的分子,不可能把所有电子均成对地排布于相应的分子轨道之中,体系中将有未配对的电子。
此时,电子体系处于开壳层状态。
显然,当电子从闭壳层分子基态占据分子轨道跃迁到基态未占据的空轨道上去,也会产生类似状态。
图2-2处理开壳层的两种方法
当描写开壳层体系时,波函数ψ一般应满足条件(2.13a)的一些Slater行列式线性组合(2.15)构成。
从而,计算方案将更加复杂。
然而,对于开壳层体系对应于极大多重度的状态来说,可以保持波函数的单行列式表示。
非限制Hartree-Fock(UHF)方法[4,9]是描写这类体系的可能方法之一。
UHF方法的基本假定是,α自旋电子所处的分子轨道不同于β电子。
从而,与闭壳层体系不同,在UHF法中引入两组分子轨道:
p个α自旋电子置于分子轨道集合ψiα中,而q个β自旋电子置于分子轨道集合ψiβ中。
电子体系的波函数为
(2.49)
由于α自旋和β自旋电子占据不同的空间轨道,所以应用ψUHF就在每种程度上考虑了不同自旋电子的相关效应。
常用的还有另一种处理开壳层体系的限制型Hartree-Fock(RHF)方法。
两种不同方案的图像表示如图2.2.所示。
把(2.49)代入(2.16)式中,可以求出ψUHF所描写的体系的能量表达式
(2.50)
其中交换积分Kijα和Kijβ分别用分子轨道ψiα和ψiβ计算。
当ψα=ψβ和p=q时,(2,50)式自动还原为(2.21)式,而ψUHF也就变为闭壳层体系的波函数(2.14)式。
如前所述,对于分子轨道ψiα和ψiβ分别向原子轨道φμ做LCAO展开有
(2.51)
此时,展开系数cα不同于cβ。
ψiα和ψiβ各自满足正交归一化条件,因为它们不同的自旋因子保证了(2.49)式中的ψiα与ψiβ的正交性。
可以分别引入α自旋和β自旋电子的原子轨道密度矩阵元
(2.52)
显然,总电子密度矩阵元等于两者的和,即
(2.53)
而两者之差定义了自旋密度矩阵元
(2.54)
将ψiα和ψiβ的展开式(2.51)代入(2.50)式中积分表达式,可得到用原子基AO表达的能量公式,并考虑到(2.52)式,E可以写为
(2.55)
相对于cμiα和cμiβ各自使用变分法极小化能量(2.55),就导出联立的两套方程组,可计算分轨道ψiα和ψiβ的能量εiα和εiβ及系数cμiα和cμiβ。
(2.56)
上式中α自旋和β自旋电子的Hartree-Fock算符矩阵元为
(2.57)
或者写成下列矩阵公式
(2.58)
正如闭壳层情况一样,方程(2.56)也是系数cα和cβ的三次联立方程组,只能用迭代法求解。
UHF方法的缺点是,波函数ψUHF满足条件(2.13b),Ms=(p-q)/2,但不是S2的本征函数,即不对应于任一个总自旋值。
在一般情况下,可以把ψUHF表成具有不同自旋多重度波函数的线性组合
(2.59)
就数值来说,展开式(2.59)中系数Cs+m很快变小,所以分离出第一个混合态S'=S+1常常足够了,可以借助于消灭算符[10]
(2.60)
用原始波函数ψUHF表达与波函数AS+1ψUHF有关的密度矩阵(2.52)、总密度矩阵(2.53)和自旋密度矩阵(2.54)式,它们是UHF方法中极重要的一些物理量,详细的表达式见引文[10]。
2.1.4开壳层体系的限制性Hartree-Fock方法
假定ψiα≡ψiβ,就可以从ψUHF得到满足(2.13a)的分子电子波函数,这就导出一个Slater波函数,它的n1个轨道均为两个自旋反平行电子占据(闭壳层),而n2个轨道为自旋相同电子单占据(开壳层)。
(2.61)
其图像如图2.2所示,波函数ψUHF是算符S2和Sz的本征函数,同时描写极大多重度状态S=MS=n2/2。
所谓限制性Hartree-Fock法就是用(2.61)型波函数作运算的。
把(2.61)代入(2.16)式,作一些变换后,可得到RHF法的电子能量表达式
(2.62)
其中k和l表示闭壳层部分的分子轨道,而m和n表示开壳层的分子轨道。
量f表示开壳层的占据程度;a和b是Roothaan常数[11],它们取决于所研究电子体系的具体特征。
例如,对半充满的开壳层体系,f=1/2,a=1,b=2。
公式(2.62)中其余量的含义同2.2节。
在LCAO近似下,此时闭壳层和开壳层分子轨道均可按原子轨道基集合φμ展开成(2.19)式。
运用关系(2.28)-(2.30)就可用AO基改写(2.62)式。
当使用矩阵表示时,可以写成
(2.63)
其中ν1=2和ν2=2f是闭壳层和开壳层的填充数;R1和R2分别是相应于闭壳层和开壳层部分的原子基的密度矩阵[见(2.33)式]。
闭壳层和开壳层部分的电子相互作用矩阵分别是G1和G2,它们是
(2.64)
矩阵G和G'可以用库仑矩阵[(2.38)式]和交换
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第一性 原理 计算 方法