第一性原理计算原理和方法

第二章 计算方法及其基本原理介绍
化学反应的本质是旧键的断裂和新建的形成，参与成键原子的电子壳层重新组合是导致生成稳定多原子化学键的明显特征。

因此阐述化学键的理论应当描写电子壳层的相互作用与重排，借助求解满足适当的Schrodinger 方程的波函数描写分子中电子分布的量子力学，为解决这一问题提供了一般的方法，然而，对于一些实际的体系，不引入一些近似，就不可能求解其Schrodinger 方程。

这些近似使一般量子力学方程简化为现代电子计算机可以求解的方程。

这些近似和关于分子波函数的方程形成计算量子化学的数学基础。

SCF-MO 方法的基本原理
分子轨道的自洽场计算方法(SCF-MO)是各种计算方法的理论基础和核心部分，因此在介绍本文计算工作所用方法之前，有必要对其关键的部分作一简要阐述。

Schrodinger 方程及一些基本近似 为了后面介绍各种具体在自洽场分子轨
道(SCF MO)方法方便，这里将主要阐明用于本文量子化学计算的一些重要的基本近似，给出SCF MO 方法的一些基本方程，并对这些方程作简略说明，因为在大量的文献和教材中对这些方程已有系统的推导和阐述[1-5]。

确定任何一个分子的可能稳定状态的电子结构和性质，在非相对论近似下，须求解定态Schrodinger 方程
''12121212122
ψψT p B A q p A p pA A pq AB B A p A A A E R Z r R Z Z M =⎥⎥⎦
⎤⎢⎢⎣⎡-++∇-∇-∑∑∑∑∑∑≠≠ （） 其中分子波函数依赖于电子和原子核的坐标，Hamilton 算符包含了电子p 的动能和电子
p 与q 的静电排斥算符，
∑∑≠+∇-=p q p pq
p
e r H 12121ˆ2 以及原子核的动能
图2-1分子体系的坐标
∑∇-=A A
A
N M H 2121ˆ 和电子与核的相互作用及核排斥能
∑∑≠+-=p A B A AB B A pA
A eN R Z Z r Z H
,21ˆ 式中Z A 和M A 是原子核A 的电荷和质量，r pq =|r p -r q |，r pA =|r p -R A |和R AB =|R A -R B |分别是电子p 和q 、核A 和电子p 及核A 和B 间的距离（均以原子单位表示之）。

上述分子坐标系如图所示。

可以用V(R,r)代表－式中所有位能项之和 ∑∑∑-+=≠≠p A pA
A B A q p pq AB B A r Z r R Z Z r R V ,1
2121),(
原子单位
上述的Schrodinger 方程和Hamilton 算符是以原子单位表示的，这样表示的优点在于简化书写型式和避免不必要的常数重复计算。

在原子单位的表示中，长度的原子单位是Bohr 半径
能量是以Hartree 为单位，它定义为相距1Bohr 的两个电子间的库仑排斥作用能 质量则以电子制单位表示之，即定义m e =1 。

Born-Oppenheimer 近似
可以把分子的Schrodinger 方程改写为如下形式
''),(21212
2ψψT p p
A A A E r R V M =⎥⎦⎤⎢⎣⎡+∇-∇-∑∑ 由于组成分子的原子核质量比电子质量大103倍～105倍，因而分子中电子运动的速度比原子核快得多，核运动平均速度比电子小千倍，从而在求解电子运动问题时允许把电子运动独立于核运动，即认为原子核的运动不影响电子状态。

这就是求解式的第一个近似，被称作Born-Oppenheimer 近似或绝热近似。

假定分子的波函数Ψ′可以确定为电子运动和核运动波函数的乘积
)(),(),(r Φr R r R ψψ=' 其中Ф(R )只与核坐标有关，代入方程有
对于通常的分子，依据Born-Oppenheimer 原理有：?A Ψ和?A 2Ψ都很小，同时M A ≈103～105，从而上述方程中的第二项和第三项可以略去，于是 易知
也即该方程可以分离变量而成为两个方程
ψψψ)()
,(212R E r R V p
p =+∇-∑ Φ
E ΦR E ΦM T A
p A
=+∇-∑
)(21
2 方程为在某种固定核位置时电子体系运动方程，而方程时核的运动方程。

E (R )固定核时体系的电子能量，但在核运动方程中它又是核运动的位能。

此时分子总能量用E T 代表。

因此，在Born-Oppenheimer 近似下，分子体系波函数为两个波函数的乘积式。

分子中电子运动波函数Ф(R )分别由和式确定。

电子能量E (R )为分子的核坐标的函数，从式看出它又是核运动的位能。

在空间画出E (R )随R 的变化关系称为位能图。

单电子近似
体系的电子与核运动分离后，计算分子的电子波函数Ψ归结为求解下面的方程
ψψE r Z r p p A pA A q p pq p
=⎥⎥⎦⎤⎢⎢⎣
⎡-+∇-∑∑∑≠'12121,2 式是量子化学的基本方程，目前已有多种求解这个方程的方法。

这些方法的区别首先是构成Ψ的方式及其相应的近似。

最常用的是Hartree 建议的单电子近似［6］。

在多电子体系中，所有电子势相互作用的，其中任意电子运动依赖于其它电子的运动。

Hartree 建议把所有电子对于每个个别电子运动的影响代换成某中有效场的作用。

于是每个电子在核电荷及其余电子有效场产生的势场中运动仅依赖于电子坐标。

从而，电子运动分开了，对于多电子体系中每个电子可以引入单电子波函数，这种单电子波函数是式单电子Schrodinger 方程的解，其中含有算符1/r pq 项，用只依赖于所研究电子坐标的有效场代替。

整个多电子体系波函数等于所有电子的单电子波函数（轨道）乘积。

电子还具有自旋角动量s ，其分量s x ,s y 和s z 满足普通角动量算符的对易关系。

算符s 2和s z 完全给定了电子的自旋，电子自旋波函数?(?)满足方程
)()(ˆξηξηz z m s
= 其中?是自旋坐标，通常把对应于自旋1/2的波函数记为?(?)，而把自旋m s =-1/2波函数记作?(?)。

在非相对论近似下和不存在外磁场时，电子的自旋和空间坐标无关，因此，因此电
子的自旋轨道可取成
)(),,(),,,(ξηψξψz y x z y x ='
考虑到自旋变量的多电子波函数由自旋轨道组成，他应当是体系总自旋S2及其Sz 的本征函数
ψψ)1(ˆ2+=S S S
ψψS
z M S =ˆ 构成体系多电子波函数Ψ时，必须考虑Ψ相对于任一对电子交换的反对称性要求，此所谓Pauli 原理［7］。

因此，一般不求出Hartree 方法的简单乘积型波函数Ψ，而是求出对应于按自旋轨道电子的所有可能置换方式的Slater 行列式波函数，此为Hartree-Fock 方法。

对于置于n=N/2轨道的Ψ上的N 电子体系，单电子近似下波函数
Ψ写为
)
()()()()()()()()
2()2()2()2()
2()2()2()2()1()1()1()1()
1()1()1()1(!
1111111N N N N N N N N N Ψn n n n n n βψαψβψαψβψαψβψαψβψαψβψαψΛΛ
Λ
Λ
ΛΛ=
该式的Slater 行列式是保证反对称性要求的唯一这类函数。

引入单电子近似便确定了波函数Ψ的形式，用它可以求解方程。

显然在一般的情况下，Ψ应当包含型行列式的线性组合，同时满足式的限制。

若式中自旋部分是单电子自旋投影算符S z 的本征值，则式就满足。

当分子的n 个轨道每个均为自旋反平行电子对占据时（闭电子壳层），一个行列式波函数就已满足和。

对于含有未配对的电子体系，这是做不到的，此时体系波函数是对应于各种轨道填充方式（不同组态）的Slater 行列式Ψl 的的线性组合
l l
l a Ψψ∑=
当适当选择行列式前系数a l 时，条件和波函数的反对称性要求均可以满足。

由于存在着电子运动的相关，不明显处理式中1/r p q 项的单电子近似，完全忽略了这种相关效应，所以，Hartree-Fock 单电子近似使波函数的计算产生了误差。

变分原理
上述单电子近似只是给出了所求解体系多电子波函数的一种形式，变分法提供了求解方程的一种方法。

Schrodinger 方程的解对应于稳定态能量。

因此若波函数Ψ是.的解，那么对于任意
微小变化?Ψ，取能量平均值
⎰>==<τΨd H Ψ|ΨH Ψ|E *ˆˆ
的变分应等于零，即
0ˆ>=<=|ΨH
Ψ|E δδ 式中积分是对Ψ的所有变量进行的，并且已假定Ψ是归一化的，即
1ˆ=⎰τΨd H
Ψ
*
由于我们寻找对应于体系基态的波函数，总能 量应当是极小值。

因此，对单电子轨道施行变分就给出这种型式波函数，能量是极小值并满足式。

从而求得的波函数Ψ就是多电子体系基态Schrodinger 方程所欲求的解。

显然，为了施行变分，波函数Ψ的型式应当充分好。

两种途径可以保证这一点：①取展开式是从充分多项，且固定轨道Ψ只对系数a l 变分；②局限于尽可能少的行列式Ψl ，若有可能做到就取一个，但此时把每个ψ表成可能的简单形式。

鉴于这种选择，区分出两类广泛应用的量子化学方法，价键(VB)法和分子轨道法(MO).
在价键法中，用孤立原子的原子轨道(AO)作为单电子波函数ψ去构成Slater 行列式Ψl 。

原子轨道的不同选择对应于不同的行列式Ψl 。

对于施行的变分，可得到确定系数a l 的方程。

为了充分靠近体系的能量，必须在式中选用足够的多项，即用多行列式波函数进行运算。

用原子轨道线性组合分子轨道(LCAO MO)法提供了另外一种选择相应于体系能量极小的多电子波函数方法。

此时，对应于分子中单电子态的分子轨道ψi 写成原子轨道φμ（基函数AO ）的线性组合
∑==m
i i c 1
μμμϕψ
实际上，这种展开有完全合理的基础。

因为靠近某个原子的电子所受的作用基本上是由该原子产生的场引起的，所以该区域中电子波函数应当近于原子轨道。

展开该式对求解变分问题的优点是明显的。

如果式中选用极大数目的项，那么VB 法和MO 法就都给出同样的能量E 和波函数ψ，当然表达式不完全相同。

这种唯一性的原因很简单，因为使用LCAO MO 的的每个行列式均可以展开为AO 组成的一些行列式。

在一般情况下，每个MO 组成的行列式应展开成AO 组成的所有行列式。

因为波函数ψ应通过AO 组成的行列式完全集合表达，从而，当使用完全集合时，MO 法与VB 法所描述的ψ就等价。

当然，不用完全及表达时，两种方法
的等价性就破坏了。

在极端性况下，某种方法中可以取一个行列式，此时可以直接看到MO 法的优越性。

对于MO 法，允许采用单行列式表达ψ（至少对于闭壳层体系），进而，通常由一些正交分子轨道组成行列式
ij j d H i
δτϕϕ=⎰ˆ*
其中δij 是Kronecker 符号。

从而是计算大为简化，并能比VB 法更简单地确定式的方程系数。

同时，MO 法的基本方程能很好的适应现代电子计算机的能力。

由于这个原因，现代的MO 方法已经成为最常用的计算多电子分子的电子结构的基本方法。

2.1.2 闭壳层体系的Hartree-Fock-Roothaan 方程
在分子轨道范围内，对闭壳层体系，在单电子近似下，用两个自旋反平行电子填充每个分子轨道ψ，可以构成一个Slater 行列式型波函数，选择轨道的自旋部分满足式，则保证了条件。

根据变分原理，若轨道ψ使得分子能量取极小值，就求出了所研究多电子体系方程的解。

将波函数代入式，并进行一些推导［见引文1－4］，可得闭壳层分子的电子能量表达式
∑∑-+=n
ij ij ij i
ii K J H E )2(2
此处H ii 是对应于分子轨道ψi 的核实Hamilton 量H core (1)的单电子矩阵元
⎰=τϕϕ(1)d ˆ(1)i
*core i ii H H 而H core 包含电子动能算符和分子中原子核对电子的吸引能算符
∑-∇-=+=A A
A core r Z V T H
12)1(21ˆˆ)1(ˆ 下面两式分别表示库仑积分J ij 和交换积分K ij
2112
**)2()1(1
)
2()1(ττϕϕϕϕd d r J j i j i ij ⎰⎰= 212**)2()1(1
)
2()1(ττϕϕϕϕd d r
K i j j i ij ⎰⎰= 积分取遍电子1和2的全部空间坐标。

从至可以看出式中各项的物理意义。

显然，单电子积分H ii 表示在核势场中分子轨道上电子能量，由于每个ψi 轨道上占据两个电子，所以乘以2。

双电子库仑积分J ij 表示
ψi 和ψj 轨道上两个电子间平均排斥作用。

由于波函数的反对称性要求出现了交换积分K ij （在Hartree 方法中不考虑它），减小了不同轨道ψi 和ψj 上平行自旋电子间相互作用，正是它们描写了相同自旋电子运动的交换相关。

然而，在Hartree-Fock 方法中，还是没有考虑反平行自旋电子间库仑排斥引起的电子相关效应。

为了求解Ψ的最优近似，必须选择一定形式的分子轨道ψi 是总能量最小。

这些分子轨道相互正交，在LCAO 近似下表成原子轨道ψμ的展开式。

Roothaan 最先解决了这一问题［8］。

关于ψi 对AO 基展开系数的方程成为Hartree-Fock-Roothaan 方程（间记为HFR 方程）。

下面简要推导这个方程。

首先从写ψi 的LCAO 展开式
分子轨道正交归一化条件给出对MO 系数的附加限制
ij
j
i
S c c δμν
μννμ=∑*
其中S μν是原子轨道ψμ和ψν间的重叠积分
⎰=τϕϕνμμνd S *
我们不难用AO 基写出至式
∑=μν
μννμH c c H i i ii *
∑><=μνλσ
σλνμλσμν|*
*j
j
i
i ij c c c c J
∑><=
μνλσ
σλνμλνμσ|**j
j
i
i ij c c c c K
其中H μν是核实Hamilton 量相对原子轨道ψμ和ψν的矩阵元。

双电子相互作用积分<μν|λσ>表达式为
21*
12
*
)2()1(1)
2()1(|ττϕϕϕϕλσμνσλνμd d r ⎰⎰>=
< 它表示电子云分布ψμψν和ψλψσ间相互作用。

从而能量E 的表达式变为
]||2[2**
*
**∑∑∑∑∑><-
><+=μνλσ
σλνμμνλσ
σλνμμνμννμλνησλσηνj
j
i
i
ij j j i i occ
i i i c c c c c c c c H c c E
引入原子轨道基的电子密度矩阵元
∑=occ
i i i c c P νμην*
2
上式化为
]|2
1
|[21><-><+
=∑∑λνμσλσμνμνλσλσ
μνμν
μνμνP P H P E 或者令
μνμνP R 2
1
=
则可以用用矩阵形式写为
RG S RH S E P P +=2
其中电子相互作用矩阵G 定义为
)()2(R K R J G -=
矩阵J (R )描写库仑相互作用，而K (R )描写电子交换相互作用，其矩阵元分别为
∑><=λσ
λσμνλσμν|)]2([R R J
∑><=λσ
λσμνλνμσ|)]([R R K
作用于两个矩阵的运算Sp 是指把此两个矩阵所响应矩阵元的乘积求和。

在LCAO 分子轨道法中，通常假定基原子轨道是固定的，而型式不变，因此对分子轨道变分归结为对展开系数c μν的变分
∑=μ
μμϕδδψ)(i i c
应当指出，式中的ψμ不变性在原则上是不必要的，已有把ψμ看成可变的一类计算方法，然而是极其困难和复杂的，常用的MO 法一般不做这类计算。

当对分子轨道变分时，借助Lagrange 乘法可以是能量极小化；这是极小化泛函 即变化系数c μi ，使稳定点δG＝0。

其中E 由方程定义。

取G 的变分，则有
}]|21
|[{2 2]| |2)[(2ˆ,
***
***=+-><-><+
=+-><->
<++=∑∑∑∑∑∑∑∑∑∑∑共轭复数共轭复数
j
j ij j i
j
j
i i
i ij
j i ij i
ij j i j i j i j i i i S c c c c H c c S c c c c c c c c c c H c c G μννν
λσ
νσλμννμμμν
μννμμνλσ
σνλμσνλμμνμννμελνμσλσμνδδελνμσλσμνδδδδ
由于δc μi *的任意性，所以必须有
∑∑∑=-><-
><+ν
λσ
μννλσμνελνμσλσμν0)]|21
|([j
j ij S c P H
因为分子轨道在酉变换下的确定性，可以自由选择非对角Lagrange 乘数为零。

定义Fock 矩阵元
)|2
1
|(∑><-
><+≡λσ
λσμνμνλνμσλσμνP H F 或记为矩阵形式
G H F +=
其中H 为Hartree-Fock 矩阵的单电子部分，而G 为双电子部分。

于是方程可以写为
0)(=-∑ν
νμνμνεi
i
c S F
或记为矩阵形式
SCE FC =
此为HFR 方程，其中E 是Hartree-Fock 算符本征值组成的对角矩阵，S 是基原子轨道的重叠积分矩阵。

解或，就给出按原子轨道展开的分子轨道系数和单电子分子轨道能量εi 。

然而，上述HFR 方程是数学上广义特征值问题。

为了求解方便，先借助于下述变换将它化成标准特征值问题。

由于S 是Hermite 矩阵，可以用酉变换θ使其对角化，即
⎥⎥⎥⎥
⎦
⎤⎢⎢⎢
⎢
⎣
⎡=+nn d d d D S O
22
11
θθ 从而可由其对角元平方根倒数构成矩阵S-1/2，同时S1/2S-1/2和S-1/2SS-1/2=1成立。

此时可做变换
2
12
1-
-
=FS
S
F T
C S C T
2
1=
则方程化成 亦即
E C C
F T T T =
这是标准特征值问题，可以采用数学上标准的对角化方法处理。

显然这些方程都是系数c μi 三次联立方程组，必须用迭代法求解：当F (0)=H 时，解上述方程或就得到MO 系数的零级近似C (0)。

用C (0)计算储F (1)；再代入方程确定新的系数C (1)，然后再计算F (2)等等。

最后，当前后两次迭代所得系数C (n )与C (n-1)符合收敛精度时，迭代
过程收敛。

应用现代快速计算机易于实现这样的计算过程。

2.1.3 开壳层体系的非限制性Hartree-Fock 方法
对于含有奇数电子的分子，不可能把所有电子均成对地排布于相应的分子轨道之中，体系中将有未配对的电子。

此时，电子体系处于开壳层状态。

显然，当电子从闭壳层分子基态占据分子轨道跃迁到基态未占据的空轨道上去，也会产生类似状态。

当描写开壳层体系时，波函数ψ一般应满足条件的一些Slater 行列式线性组合构成。

从而，计算方案将更加复杂。

然而，对于开壳层体系对应于极大多重度的状态来说，可以保持波函数的单行列式表示。

非限制Hartree-Fock(UHF)方法[4,9]是描写这类体系的可能方法之一。

UHF 方法的基本假定是，α自旋电子所处
的分子轨道不同于β电子。

从而，与闭壳层体系不同，在UHF 法中引入两组分子轨道：
p 个α自旋电子置于分子轨道集合ψi α中，而q 个β自旋电子置于分子轨道集合ψi β中。

电子体系的波函数为
)}
()()12()12()2()2( )3()3()2()2()1()1(det{)!
(11q 211q p q p q q q q q p Ψq p q UHF +++++=
++αϕαϕβϕαϕβϕαϕα
αβα
βαΛΛ
由于α自旋和β自旋电子占据不同的空间轨道，所以应用ψUHF 就在每种程度上考虑了不同自旋电子的相关效应。

常用的还有另一种处理开壳层体系的限制型Hartree-Fock (RHF)方法。

两种不同方案的图像表示如图.所示。

把代入式中，可以求出ψUHF 所描写的体系的能量表达式
∑∑∑∑++--+=q
p i q p ij ij
ij p ij ij ij ii K K J H E )(21β
β
α
其中交换积分K ij α和K ij β分别用分子轨道ψi α和ψi β计算。

当ψα＝ψβ和p=q 时，(2,50)式自动还原为式，而ψUHF 也就变为闭壳层体系的波函数式。

如前所述，对于分子轨道ψi α和ψi β分别向原子轨道φμ做LCAO 展开有
⎪
⎭
⎪⎬⎫====∑∑ββμ
μβμβααμ
μα
μαϕψϕψϕψϕψc c c c i i i i or
or
图2-2 处理开壳层的两种方法
此时，展开系数c α不同于c β。

ψi α和ψi β各自满足正交归一化条件，因为它们不同的自旋因子保证了式中的ψi α与ψi β的正交性。

可以分别引入α自旋和β自旋电子的原子轨道密度矩阵元
∑∑==q i
i
i p
i i
i c c R c c R
βνβηβμνα
ναηαμν
**
显然，总电子密度矩阵元等于两者的和，即
βμναμνηνR R P +=
而两者之差定义了自旋密度矩阵元
β
μν
αμνμνρR R spin -= 将ψi α和ψi β的展开式代入式中积分表达式，可得到用原子基AO 表达的能量公式，并考虑到式，E 可以写为
]
|)( |[21
><+-><+
=∑∑λνμσλσμνφλσβμνα
λσ
αμνμνλσλσμνμν
μνμνR R R R P P H P E 相对于c μi α和c μi β各自使用变分法极小化能量，就导出联立的两套方程组，可计算分轨道ψi α和ψi β的能量εi α和εi β及系数c μi α和c μi β。

)(0
)(=-=-∑∑ν
βνμνβ
β
μνν
ανμνααμνεεi
i i
i
c S F c S F
上式中α自旋和β自旋电子的Hartree-Fock 算符矩阵元为
∑∑><-><+=><-><+=λσ
βλσ
λσμνβμν
λσ
αλσλσμναμνλνμσλσμνλνμσλσμν]
||[]
||[R
P H F
R P H F
或者写成下列矩阵公式
)
()()()()()(β
β
α
β
β
βααααR J R K R J H G H F R J R K R J H G H F +-+=+=+-+=+=
正如闭壳层情况一样，方程也是系数c α和c β的三次联立方程组，只能用迭代法求解。

UHF 方法的缺点是，波函数ψUHF 满足条件，M s =(p-q )/2，但不是S 2的本征函数，即不对应于任一个总自旋值。

在一般情况下，可以把ψUHF 表成具有不同自旋多重度波函数的
线性组合
∑=++=q
m m s m s UHF ΨC Ψ0
就数值来说，展开式中系数C s+m 很快变小，所以分离出第一个混合态S'=S+1常常足够了，可以借助于消灭算符［10］
)2)(1(ˆˆ21
++-=+S S S A S 用原始波函数ψUHF 表达与波函数A S+1ψUHF 有关的密度矩阵、总密度矩阵和自旋密度矩阵式，它们是UHF 方法中极重要的一些物理量，详细的表达式见引文［10］。

2.1.4 开壳层体系的限制性Hartree-Fock 方法
假定ψi α≡ψi β，就可以从ψUHF 得到满足的分子电子波函数，这就导出一个Slater 波函数，它的n 1个轨道均为两个自旋反平行电子占据（闭壳层），而n 2个轨道为自旋相同电子单占据（开壳层）。

)}
2()2()2()2( )
12()12()2()2()1()1(det{)!
2(12121211111111112n n n n n n n n n n Ψn n n n RHF ++--+=
+αϕβϕαϕβϕαϕΛΛ
其图像如图所示，波函数ψUHF 是算符S 2和S z 的本征函数，同时描写极大多重度状态
S=M S =n 2/2。

所谓限制性Hartree-Fock 法就是用型波函数作运算的。

把代入式，作一些变换后，可得到RHF 法的电子能量表达式
∑∑∑∑∑-+-++-+=mn
km
km km mn mn kl
m
mm
kl kl k
kk K J bK aJ f H f K J H E )
2(2)2( 2[)2(2
其中k 和l 表示闭壳层部分的分子轨道，而m 和n 表示开壳层的分子轨道。

量f 表示开壳层的占据程度；a 和b 是Roothaan 常数[11]，它们取决于所研究电子体系的具体特征。

例如，对半充满的开壳层体系，f=1/2，a=1，b=2。

公式中其余量的含义同节。

在LCAO 近似下，此时闭壳层和开壳层分子轨道均可按原子轨道基集合φμ展开成式。

运用关系-就可用AO 基改写式。

当使用矩阵表示时，可以写成
)2
1
()21(222111G H R S G H R S E p p +++=νν
其中ν1=2和ν2=2f 是闭壳层和开壳层的填充数；R 1和R 2分别是相应于闭壳层和开壳层部分的原子基的密度矩阵［见式］。

闭壳层和开壳层部分的电子相互作用矩阵分别是G 1和G 2，它们是
)
()()
()(2211122111R G R G G R G R G G νννν'+=+=
矩阵G 和G'可以用库仑矩阵[式]和交换作用矩阵[]表示为
)(2
1
)()(R K R J R G -
= )(2
1
)()(R bK R aJ R G -='
矩阵描述闭壳层电子相互作用，而矩阵描述开壳层电子相互作用。

用极小化式，可以得到确定分子轨道展开成AO 基的系数方程，当然这必须在分子轨道相互正交的条件下进行，分别对闭壳层和开壳层轨道的LCAO 系数变分。

通常这导出两个型的Hartree-Fock 矩阵，对闭壳层和开壳层分别有
2
22111G H F G H F +=+=
还可以由不同壳层分子轨道正交条件导出两个型的方程组。

然而，正如Roothaan 指出的［11］，对于所研究的函数，可以将两个型的方程组统一成一个求本征值和本征向量的方程
SCE FC =
当然，这需要按下述公式确定上式中的Hartree-Fock 矩阵
322113121212
)(R F F R R F R R F R F '-'+''+''=νν 矩阵R i '=1-R i (i =1,2)，R 3'=1-R 1-R 2，此处1为单位对角矩阵。

可以把式写成与式相应的形式
B BR RB K J K J H F -++-+-+=221122
这里B =2aJ 2-bK 2；R 是总密度矩阵R =R 1+fR 2。

库仑积分和交换积分下角标1和2分别对应于闭壳层和开壳层，包含有式定义的Hartree-Fock 算符方程的求解归结为逐次迭代法。

矩阵R 1由n 1个闭壳层轨道计算，而R 2由n 2个开壳层轨道计算。

波函数式以及公式-不仅可用于描写极大自旋多重度状态，而且也可用于描写开壳层轨道部分填充的一些状态。

此时，只须改变式中a ，b ，f 数值，也并不使问题复杂化。

Roothaan 论文中已给出用RHF 方法描述分子状态的这类系数［11］。

Pitaer 等人推广了Roothaan 的研究，除了给出原子和线型分子能量系数外，还给出了Td 和Oh 分子的能量系数，使该方程更一般化了［12］。

2.1.5 SCF 计算结果的分析
通过SCF 计算得到了收敛波函数|ψ>，就可以应用此波函数|ψ>分析计算常用的物理量，现将主要的物理量计算介绍如下。

总能量
体系的电子总能量E 是Hamilton 算符的平均值<ψ|H |ψ>，即
∑∑∑+=
-+>='=<2
/,2
/)(21)2(2||N j
i core
ij ij N i
ii F H P K J H ΨH ΨE μν
μνμνμν 而体系的总能量E T 是电子能量加上原子核－原子核排斥能量
∑∑
>+=A
A B AB
B
A T R Z Z E E
这是最有意义的量，因为预期稳定分子的平衡几何构型在E T 的极小点处。

Koopmans 定律
电离势的理论计算式是提供说明分子光谱的有用工具。

与半经验计算不同，从头算处理了全部电子，因此不必区分价电子和内壳层电子电离势。

一般不去区分光电子光谱（PES ）和学分析电子光谱(ESCA)。

两者在理论上是统一处理的。

Koopmans 定律提供了电离势的理论上最简单近似[13]，它使亲态闭壳层体系的Hartree-Fock 轨道能量εh 的负值等于第h 个电离势
h K h IP ε-=
所以，按Koopmans 定律，若理论上有若干个电子占据分子轨道，则光电子光谱应当存在有相应个谱带，粗略地说，在多数情况下实验观察的确如此[14]。

然而，定量地讲，方程并非很满意地列数了观测谱。

这是由于Koopmans 定律有三个内在缺陷：①Koopmans 定律不管从分子到离子时轨道的自洽重排，亦即忽略了与这种重排有关的能量ΔR ，通常叫做重新组合能或弛豫能。

总之Koopmans 定律认为从分子到离子所有占据的分子轨道保持不变，这是一种近似。

②没有考虑从分子到离子的相关能变化ΔC 。

③假定分子和离子的的相对论能量相同。

首先，相对论能量随着原子系数的增加而增大，并且对于内壳层来说是很大的数值。

例如对内壳层电离势来说，CH 4的相对论修正为，对于Ne 的修正为，而对于Ar 对修正则是14eV 。

由于从头计算不直接计算相对论修正。

其次，当忽略相对论效应时，对于内壳层从头计算给出电离势的估计值数值偏低。

从而不考虑相对论效应时，电离势应为
C R IP h KRC h ∆+∆+-=ε
差值IP h KRC -IP h K 是所谓的Koopmans 缺欠，即应当考虑电子弛豫和相关能修正。

当采用ΔSCF 方法时，可以考虑电子重排能ΔR
SCF SCF h KR h E E IP -=
其中E SCF 是电离前分子的SCF 总能量，而E h SCF 是空穴态离子的SCF 能量。

这表明，电离势是绝热SCF 计算的分子和离子的总能量之差，有时称之为竖直电离能。

单电子算符的平均值
诸如分子的偶极矩、四极矩、核位置的场梯度、反磁化率等可以采用单粒子算符的和描写，一般形式为
∑==N
k k g
O 1)(ˆˆ 其中g (k )是任意一种只依赖于单电子坐标的算符。

单粒子算符的平均值是
∑∑><>=<>=>=<<μν
νμμνϕϕϕϕ
ψψ|ˆ||ˆ||ˆ|ˆ2/11
1g P O O O N i
i
i
从而，除密度矩阵外，只须计算单电子积分<ν|g |μ>。

下面以偶极矩为例加以说明。

偶极矩的经典定义为
∑=i
i i r q μ
向量r i 除的电荷为q i 。

偶极矩相应的量子力学定义为
∑∑+>-=<A
A A N
i
i R Z r ψψμ||
其中第一项是电子（电荷为-1）对偶极矩的量子力学贡献，第二项是原子核的核电荷Z A 的经典贡献项。

电子偶极矩算符∑-N
i i r 是单粒子算符。

所以应当有
∑∑+><-=A
A A R Z r P μνμν
μνϕϕμ||
这是一个向量表达式，其分量公式可以写为
∑∑+><-=A
A A x X Z x P μνμν
μνϕϕμ||
对于μy 和μz 也有类似的公式。

为了计算偶极矩，除了密度矩阵外P μν外只须计算偶极矩积分
⎰>=<1111*)()(||dr r x r x μνμνϕϕϕϕ
至于与电荷密度有关的问题下面加以说明。

集居数分析
在分子轨道计算中，人们通常把分子中电荷密度定义为
∑=μν
νμμνϕϕρ)()()(11*
r r P r i
它表示空间任意点发现电子的几率，常常画出通过分子中各种截面的电荷密度轮廓图[16]。

但是，由于分子给定原子或原子核中的电子数不是唯一确定的，有时可以给出不同的集居数分析。

因为
∑⎰=2
/21|)(|2N i i r dr N ψ
每个MO 占有两个电子。

把ψi (r )的展开式代入上式就有
PS
T PS S P S P r r dr c c r r c c dr N r i
i i N i i i ======∑∑∑∑∑∑∑∑⎰∑⎰∑∑μ
μμμ
ν
νμμσμ
ν
νμ
μσμ
ν
νμνμμν
νμνμϕϕϕϕ)( )
( )
()(2 )
()(2**2
/**
其中(PS )μμ为联系于ψμ轨道的电子数，这称为Mulliken 集居数分析[17]。

对于给定原子A 在所有轨道上的电子数为∑μ
μμ)(PS A ，从而A 原子的净电荷数为
μμμ
)(PS Z Q A A A ∑-=
其中Z A 是A 原子核的核电荷数，求和取遍所有属于A 原子的轨道ψμ。

但是，上述定义不是唯一确定的，因为T r AB=T r BA ，就有
∑-=μ
μμαα)(1PS S N
对于任意a 成立。

当a =1/2时，得到
∑∑'==μ
μμ
μ
μμP PS S N )(2
121
其中P'μμ是对称正交化基组表象中的密度矩阵
∑∑'''=μ
ν
νμμν
ϕϕρ)()()(*r r P r。