计算物理方法(Sec1, Sec2)

1、计算物理熊诗杰（用于凝聚态物理） 凝聚态物理中数值计算的基本特点凝聚态物理中数值计算的基本特点 大型矩阵的数值严格对角化方法大型矩阵的数值严格对角化方法 密度泛函理论和固体能带理论密度泛函理论和固体能带理论 Monte Carlo方法方法第一章 凝聚态物理中的数值计算凝聚态物理中的数值计算 一、凝聚态系统 基本特点 1. 由大量粒子组成，是一个多粒子体系； 2. 这些粒子都是微观粒子，满足量子力学的运动规律。 因此，凝聚态系统的状态要用多体波函数表达。 二、多体波函数满足的基本方程 - 薛定谔方程 这里， 是多体波函数，其自变量是时间和描写各个粒子状态的广义坐标， 是系统的哈密顿算符。 如果

2、不依赖时间，有我们有如下的定态薛定谔方程： , iHtHH,Eti,EHE 是能量的本征值。 三、薛定谔方程的基本特点 1. 是关于未知波函数的线性齐次方程 -如果选取一组正交归一完备的基函数，可以将波函数表达为向量，哈密顿算符表达为矩阵。 2. 只有 取一定的值才有非零的波函数 - 能量本征值 E 3. 是厄米算符，能量本征值为实数。 4. 如果哈密顿表达为 阶矩阵，则有 个能量本征值，对应着 个正交归一的本征波函数。本征值有可能在数值上一样（简并）。NN NHNHE 求能量的本征值和对应的本征波函数。波函数可以是单体的波函数，也可以是多体的波函数。数值计算的基本出发点，是构造一组完备（或近

3、于完备）的基函数，使波函数成为向量，哈密顿成为矩阵。四、数值计算的基本任务第二章：大型第二章：大型Hamiltonian系统的系统的数值严格对角化方法数值严格对角化方法Lanczos方法方法第一节 概述 在对一个实际的物理系统进行理论研究的过程中，经常会遇到要求一个大型或超大型Hermit矩阵的本征值和本征函数的问题。矩阵的维数往往会达到几十万甚至上百万。这样的矩阵在包含有大量微观自由度的宏观体系中是经常遇到的。 例如，如果我们用分立的自旋来描写磁性晶格中的磁性原子，而磁性原子间的互作用则用相邻的自旋格点间的交换积分来表达，我们可以得到如下的Heisenberg哈密顿量： ,()xxyyzzx

4、yzijijiji jHJJJS SS SS S这里Jx ,Jy , Jz是三个方向上的耦合常数，Sxi , Syi , Szi 分别是第i个格点上的自旋在x, y, z三个方向上的分量。如果我们要将这个哈密顿量用数值的方法对角化，以求出其本征值和本征函数，我们就必须首先选择Hilbert空间中一组完全的向量 集作为基，把哈密顿量在这个基上用矩阵的形式表达出来，再对该矩阵进行数值求解。然而，与能带论中所用的单电子近似不同的是，这种强关联的多粒子体系所对应的Hilbert空间的基矢的个数是极其庞大的。在我们这个例子中，每个格点上的自旋的状态数是(2S+1)，如果我们研究的有限的系统包含有N个格点

5、，那么总的状态数就是(2S+1)N个，也就是说，哈密顿矩阵的大小应该是(2S+1)N(2S+1)N。这是一个极其庞大的数字，例如，我们研究一个只 有44大小的自旋为1/2的二维晶格，那么哈密顿矩阵的大小就是6553665536，这样大小的矩阵的直接数值求解在现在的绝大多数计算机中都是办不到的。 矩阵维数的增大使得求解本征值和本征函数变得非常困难，但是，我们看到，在绝大多数在物理上感到兴趣的问题中，哈密顿矩阵都有如下特点： (1) 哈密顿矩阵都是厄米矩阵。 (2) 尽管哈密顿矩阵异常庞大，但是其中绝大多数非对角矩阵元都是零，只有数量很少的非零矩阵元，也就是说，这种矩阵是极为稀疏的矩阵。例如，对我

6、们刚才讲的例子，哈密顿矩阵的每一行尽管有65536个元素，但其中非零的矩阵元不会超过5个。 (3) 所考虑的物理系统，一般都具有一定的对称性，满足一定的守恒律。这样，我们就可以根据对称性和守恒律的要求，采用一定的基函数，使哈密顿矩阵分解成为互不关联的子矩阵，从而可以对比较小的子矩阵分别求解。 Lanczos方法就是根据上述特点所发展起来的一种极为有效的哈密顿矩阵求解本征值和本征函数的数值方法。近年来，它已被广泛应用到强关联系统、无序系统等极为重要的物理系统之中，成为计算物理中一项至关重要的新方法。我们这一单元，就是介绍这一方法的基本做法和它的一些重要的应用。第二节 基本方法 发展Lanczos

7、方法的最重要的目的，就是要找出一条有效的途径，来解决哈密顿矩阵过于庞大的问题。问题的关键就在于，如何充分利用前面说到的三个特点，即厄米性，稀疏性和对称性，来达到我们的目的。这是一个极其复杂的问题，涉及到一个物理系统的一些非常具体的特点。在这一节中，我们将着重介绍如何利用体系的厄米性来将一个庞大的哈密顿矩阵约化为较小的矩阵，这是对各种物理系统都通用的基本的Lanczos方法。至于利用稀疏性和对称性来作进一步的约化，我们将在以下几节中加以介绍。 我们设H是某一在N维Hilbert空间中的厄米算符，而v1是该空间中的一个任意的归一化的向量，即：v1+v1=1。我们连续地用H去作用，来构造出一组正交归

8、一的向量：vvvvvvvvvvvHHH4333223322211221111. 我们注意到这些方程的结构是三对角型的，即每个方程的右边只出现相邻的三个元素。这是由哈密顿矩阵的厄米性所决定的。这一迭代序列将在正交归一向量走遍整个Hilbert空间后自动中止，也就是说，在做到第N步时，我们有： vvvvNNNNNNNH111但由于在Hilbert空间中最多只能有N个相互正交归一的向量，因此新的向量vN+1必须为零。这样，这一方法具有一种十分有用的性质，即：即使我们一开始并不知道向量空间的维数，这一方法也会很好地自动中止。注意到我们这里处理的是一般的算符，而不是矩阵，因此空间的维数不需要明显地写出来

9、。 这些向量，v1, ,vN , 称为Lanczos向量，它们构成了Hilbert空间的一组正交归一的完备集，以它作为基，哈密顿算符可表达为三对角矩阵的形式：.33222111HH 如果U是H的本征向量，而u是 的对应的本征向量，它们之间的关系是： UVu 这里V是以Lanczos向量作为列的矩阵。 如果算符H具有简并的本征值，则初始向量v1将包含对应于这个简并的本征值的一个本征向量（这是因为v1在这个简并本征值所对应的子空间的投影只有一个线性无关的向量）。这意味着Lanczos过程将在小于N步的地方中止，而且这个过程只给出对应于这个简并集合中的一个本征值和一个本征向量。一般来讲，如果H具有n

10、个不同的本征值 H（nN）E1，E2，, En，它们分别是p1, p2, ,pn重简并，那么上述的迭代过程将中止在第n步，所得到的由Lanczos向量所张开的n维子空间将只包含每个简并集合中的一个本征态。所得到的nn阶三对角矩阵将具有n个不同的本征值。于是，我们需要再去求出剩下的N-n个本征值和本征向量。我们可以选择一个新的初始向量，使它与所得到的所有Lanczos向量正交，来继续上述过程。同样，这个过程将经过m步之后中止，m应等于所剩下的具有不同数值的本征值的个数。我们可以这样继续下去，直到得到整个Hilbert空间中的全部本征值和本征向量为止。这样得到的三对角矩阵具有块对角形的结构，即：

11、这里每个块都是非简并的，而且其维数等于在前面的块构成之后所剩下的不相等的本征值的个数。 xxxxxxxxxxxxxxxxxxH Lanczos方法具有一个十分有趣而且至关重要的性质，即：在经过p - 1步后，所得到的p个Lanczos向量可以完全决定三对角矩阵的左上方的大小为pp的部分的元素，如果我们对这个pp的矩阵求本征值和本征向量，那么在p增大的时候，这些本征值和本征向量将会极快地收敛于完全的矩阵的低端的本征值和本征向量。对一个物理系统，在多数情况下，我们只对能量较低的本征值和本征向量感到兴趣，因为它们对应于系统的基态和能量较低的激发态。这样，我们就可以根据Lanczos方法的这个性质，不

12、去管原先的矩阵H有多大，只要做有限的p步，得到pp大小的矩阵，求出它的本征值和本征向量，然后增大p的值，使得到的低端的本征值收敛。只要它们 收敛，我们就可以认定所得到的低端本征值和本征向量就是完全矩阵的低端本征值和本征向量，从而得到所求系统的基态和较低的激发态。 上述的Lanczos方法的这一性质，已经得到了严格的证明。在计算实践中，收敛速度之快令人吃惊。在许多情况下，不管的大小有多大。所取的Lanczos步数往往只要取50到100之间，收敛的精度就可以达到计算机的精度。但这并不是说Lanczos方法很容易做。实际上，Lanczos方法是计算物理中技巧性非常高的工作。这一方面是因为所研究的物理

13、对象极为复杂，有时光是为了写出矩阵元就要费很大的工夫。另一方面，在算法中某些不引人注目的误差，如舍入误差等，都会破坏整个程序的运行。因此，上述的Lanczos方法只有在运用精确的算法和程序时才是有意义的。在实际的计算中，计算误差将以各种方式 影响到整个过程。这里我们暂不讨论Lanczos方法本身的精确性和稳定性，但我们指出一些不精确的算法所产生的后果，初看起来，这些后果似乎是十分奇怪的。例如，由于计算过程中产生的舍入误差，使得在N1步后，所得到的vN+1并不准确为零，这样整个过程就不会在N1步后中止。尽管这里的vN+1并不大，但已足以妨碍整个过程的正常运行。如果继续进行下去，同样的本征值将会被

14、一次又一次地产生。类似地，虽然从理论上来说在第N1步中vN必须与前面所得到的所有向量正交，然而在实践上vN是通过以下的公式计算的： 在本式中由于可能出现两个相近的数字相减，会引起较大的有效数字的丧失，最后的结果可能是非常奇异的。标准的方法是：为了保证其正交性必须将当前的向量v与前面的所有的向量从新进行正交化。即，必须在程序中引入如下的步骤：vvvvwnnnnnnnnH112211 这是重正交化手续，它能保证Lanczos方法能正常中止。 Lanczos方法使得我们能够把一个很大的矩阵转化为一个等效的、维数较小的三对角矩阵。为了求出对应的本征值和本征态，我们需要一个xxxvwvbvbwxnnnn

15、niiiininn211, 有效的数值方法，以便对三对角矩阵进行对角化。在下一节，我们将给出两种较为常用的方法，对三对角矩阵进行对角化。第三节 三对角矩阵的对角化 在这一节中，我们研究两种对三对角矩阵进行对角化的方法。下面我们将分别对它们进行阐述。 (1) Dean方法。 在Lanczos方法中，我们得到了有限维数的三对角矩阵，由于初始的Hamilton矩阵是Hermit矩阵，所以得到的三对角矩阵也是Hermit矩阵。这种三对角的Hermit矩阵可以很方便地用Dean方法来进行对角化。这种方法是基于如下的负本征值定理：定理如果矩阵是M阶的三对角Hermit矩阵，则它的比某一实数 小的本征值的个

16、数等于所有的实数U(i) (i=1, 2, 3, ., M) 中的负数的个数，这里实数U (i) 由如下的递推关系来决定：) 1 () 1 (,.,4 , 3 , 2),1(| ) 1(|)()(2EUMiiUitiEiUE(i) 是矩阵 的第i个对角矩阵元，而t(i) 是矩阵 的第i个次对角矩阵元。 HH 根据这个定理，我们可以方便地通过改变 的值求出我们所需要的能量范围里本征值的分布情况，并且可以求出各个本征值的准确位置。在求出本征值后，我们可以通过如下的方法求出所对应的本征函数： 我们设定对应于本征值 的本征函数在Lanczos向量所张开的空间中的第n个分量是1，则其它的分量可以通过如下

17、的递推关系求出：je) 1 (1,)(1),11(| ) 1(|)(1),1(,)(112) 1(jjMijiininineEeMEiMiteiEinMainta 最后再通过归一化，即可得到对应的本征函数。实际上，一旦我们原来的厄米矩阵已通过Lanczos方法约化为三对角形，有关确定其本征值的可能的方式就是直接寻找特征多项式pN()的根。上述方法就是通过递推过程产生出低次多项式的Sturm序列，该序列将本征值局限于一些实轴区间上，然后可以运用象二分法或Newton法那样的求根方法来精炼区间。对应的本征向量于是也可以用上述的逆迭代法求得。 （2） 运用标准的Lapack或Linpack子程序计算

18、三对角矩阵的本征值和本征向量。参考文献：E. Dagotto, Rev. Mod. Phys. 66, 763 (1994)习题：对二维Hubbard模型，用Lanczos方法求解基态能量和基态波函数：0,ijiiiii jiHt a aUa a a a这里，晶格是2x2的方型晶格，参数01,6.tU提示：1. 考虑对称性：总粒子数N守恒，总自旋z方向分量守恒。 和 分别守恒 给定总粒子数 ，则有如下子空间 ，1个基函数 ，16个基函数 ， 36个基函数 ， 16个基函数 ， 1个基函数选择基函数：对 子空间，二个自旋向上的电子有如下6个组态： 同样，二个自旋向下的粒子也有上述6个组态该子空间

19、的36个基函数可由以上两种组态直乘而得。 基函数的标号：可由两个指标 标号，每个指标可取1，2，6个值，分别为两种自旋组态的标号。NN4N 0,4NN1,3NN2,2NN3,1NN4,0NN2,2NN12 , 13 , 14 , 23 , 24 , 34, i j编制子程序进行Lanczos计算Hx第二部分 能带理论熊诗杰 第一节 概述 由于在固体中传导电子的性质决定了绝大部分实验可观察的效应的基本特征，因此有必要知道固体的电子能带结构，也就是在不同的能带中电子的能量和波矢的依赖关系，这里是能带的标号，是波矢。现在我们的能带论知识已经使我们能够有效地预言和解释大部分半导体和金属的性质，而且在新

20、型人工技术材料和新器件的研究方面发挥着重要的作用。 能带论的基础是Born-Oppenheimer 近似和单电子近似。在此近似下晶体中电子的能谱可由薛定谔方程决定：)()()(2202rrrUm这里 是晶体中作用在电子上的周期势。于是，求解方程（1.1）的问题归结为两个方面：（1） 决定周期势 ；（2） 在此势下对（1.1）求本征值和本征函数。然而，晶体中单个电子所感受到的周期势是由离子和其它电子所产生的，因此 一般需用自洽的Hartree-Fock方法来决定。尽管对此仍有一些争议，但现在用来决定周期势的大部分方法依然是在Hartree-Fock的基础上加某些近似和改进。 )(rU（1.1）)

21、(rU)(rU 如果在没有互作用的情况下单电子的波函数为 ，这里 x 包含了电子的空间坐标和自旋，则在考虑了互作用的情况下 n 个电子系统的Hartree-Fock 方程可写成：)(xi 11112222*21112222*211)()()()()()(xxrdxxxexrdxxxexHiijnjijinjjji(1.2)这里， 221010(/2)( ),HmUr )(0rU是离子产生的势， |2112rrr。等式左边的第二项是库仑项，而第三项是交换项。 Hartree-Fock方程与单电子的方程有着本质的不同：在交换项，能量所乘的是 ，而不是 ，说明对交换项，电子不能作用于其本身。jiHa

22、rtree-Fock方程可以写成等效的形式： 11111*122212*1*21112222*211)()()()()()()()()()(xxxxrdxxxxxexrdxxxexHiiinjiiijjiinjjji这个方程类似于单电子的形式。交换势不仅与坐标 有关，也与函数 有关。我们可定义交换电荷密度：1xinjiiijjixxxxxxe111*212*1*)()()()()()(如果我们把交换电荷密度对 积分，并计及 的正交性，可得到 。如果 逼近 ，交换电荷密度变成2xie2x1xnjjjxxe111*)()(以上两个结论对所有的波函数都成立。这使得我们可以取一些与波函数无关的近似的交

23、换势，使得薛定谔方程更容易求解。其中一种近似是取加权平均的交换势。权重可取为在 上的电子处于态 的几率：1xi)()()()(11*11*xxxxkkkii应用这一权重，平均交换电荷密度变为：nkkknjkkjjkxxxxxxe111*1,212*1*)()()()()()(而Hartree-Fock 方程可写成单电子的形式： 11111*1,122212*1*21112222*211)()()()()()()()()()(xxxxrdxxxxxexrdxxxexHiiinkkknjkkjjkinjjji 另一种由Slater提出的近似方法更为简洁和常用。它假设平均交换势可用同样电荷密度的自由

24、电子气的交换势来代替： ,/,11ln4121)(),(6max2023/2ppFFreUex这里 是电子的最大动量， 是单个电子所占据的半径：maxp0r3/10|43er对所有的占据动量态求平均， 成为3/4，于是：)(F3/102023/2833643nereUex这里 是电子密度。Slater方法的关键在于用局域电子密度取代 ，即0n0n3/128)(33xneUex于是，薛定谔方程变为： 113/1111*21112222*211)()()(833)()()(xxxxexrdxxxexHiiinjjjinjjji在能带计算中，一般采用原子单位制，薛定谔方程可写成： xxxxerrdx

25、xxerrZiiijjjjtNjjjj)()()(836| |) () (2|23/11*211*22maxmax这里 是第 个原胞的第 个原子的坐标， ， 是系统中的原胞数， 是原胞的体积。 是原胞中原子的个数， 是第 个核的电荷数。 是电子数。我们假设一半电子自旋向上。电子自旋数的小的偏离只会造成交换能的小的修正。 rnrn0NN0tZtZNj1max在原子单位制，能量单位是Rydberg 20113.6049 eV2eRyda12 ,21 ,12020emem第二节 密度泛函理论一、 密度泛函理论的基础HohenbergKohn定理 （一）定理一 对于不计自旋的全同电子系统，其哈密顿写为

26、 VUTH（1）其中，)()(rrdrT)() () ()(| |121rrrrrrdrdrU( )( )( )Vdrv rrr这里 和 分别表示在 r 处产生和湮灭一个粒子的费密子场算符。则其基态能量是粒子数密度函数 的唯一泛函。( ) r( ) r| )()(|)(rrr证明： 首先证明：不可能存在另外一个 ，也具有同样的密度函数 。 设含有 的哈密顿 的基态为 ，基态能为 ；含有 的哈密顿 的基态为 ，基态能为 。于是，根据变分原理，有)( rv)()( rr)(rvE)( rvHEH)()( )(| | | rvrvrdrEVVHHHE同样，我们有 )( )()(| | | |rvrv

27、rdrEVVHHHE综合上式，有 ，这是不可能的，因此， 是 的唯一泛函，进而 是 的唯一泛函，所有其它性质，包括基态能量、波函数等，都是 的唯一泛函。EEEE)(rv)(rH)(r)(r (二) 定理二 基态能量泛函 在粒子数不变的条件下对正确的粒子数密度取极小值，该极小值即为基态能量。 证明：对于全同电子系统的哈密顿（1），能量泛函 的定义为：)(rE)(rE|)()(UTrrdrvE再定义另一泛函为|UTF它不包含外场的贡献。根据变分原理，粒子数不变时，任意态 的能量| | 0VUTE在 为基态 时取极小值。设 是与外场 相联系的系统基态，则根据定理一， 和 均依赖于系统的密度函数 ，于

28、是 是 的泛函。根据变分原理，有： )( rv)( rv)( r 0E)( r)()()( )( | | 00GGErrdrvFErrdrvFEVUTE这样， 也小于所有其它与外场 相联系的密度函数 所对应的能量。亦即能量泛函的极小值为基态能量。GE)( rv)( r在没有外场时，能量泛函 可分成如下几项 |ETUV 01( ) ( ) 2|xcrrETUdrdrErr其中第一项为动能，第二项为原子核的势能，第三项为电子之间的库仑能，第四项为交换关联能。 二、 单电子化：KohnSham方程 根据上述定理，基态能量和基态的密度函数可由能量泛函对密度函数求变分得到。通过变分，得到如下方程0)()

29、(| |) ()()()()(0rrErrrdrrurrTrdrxc加上粒子数不变的条件， ，有 ( )0drr)()(| |) ()()()(0rrErrrdrrurrTxc这里拉格朗日乘子 具有化学势的意义。上式左边的后三项具有有效势的形式。但对多粒子系统的动能仍然一无所知。为此，Kohn和Sham提出，可将多粒子系统的动能泛函 用无互作用系统的动能泛函 代替，该无互作用系统具有与多粒子系统同样的密度函数，而把 和 的差别归入 ，sTTsTxcE这样就把所有复杂的因素都归入 ， 则再通过一定的假设给出。对于无互作用系统，其粒子数密度为xcExcENiirr12| )(|)(其动能泛函为Ni

30、iisrrdrT12*)()(对 的变分可用对 的变分代替，拉格朗日乘子用 代替，有)(riiE0)(/1)()()(1*rrrdrErEiNiii于是可得 2 ( ) ( )( )KSiiiVrrEr 这里，)()(| |) ()()()()()(00rrErrrdrrurVrVrurVxcxcCoulKS这就是单粒子化的KohnSham方程。 三、 局域密度泛函近似交换关联势的近似形式 交换关联势 一般是非局域的，具体计算中，KohnSham提出局域密度泛函近似。xcE（一） 交换势的近似形式1. Slater平均交换势近似3/128)(33reVex这是借助于电子气的交换势的形式，并对整

31、个Fermi球求平均得到。2. Kohn-Sham-Gaspa交换近似3/128)(32reVex这是借助于电子气的交换势的形式，但对Fermi面求平均得到。3. Slater交换能近似3/128)(33reVex13/2这是介于Slater平均交换势近似和Kohn-Sham-Gaspa交换势近似之间的形式。（二）关联势的近似形式1. Wigner关联势近似 2)79. 7(79. 73/288. 0sscrrV31) 3/4(sr2. Hedin-Lundqvist关联势近似scrV211ln2045. 0 （三） 局域密度泛函的局限性 1 不能用于强关联问题 2 对激发态的结果不好。 3

32、对d电子、f电子等结果不好。 四、 Bloch定理在能带计算中的应用 对于完整晶体，其势能具有晶格的周期性，即在KS方程)()()(2rErrViiiKS中，如果在局域密度泛函近似下把 看做 的简单函数 ，有)(rVKSr)(rVKS)()(rVRrVKSlKS这里 是晶格的格矢。对这样的周期势，根据Bloch定理，其波函数可写成lR),(),(rkuerknikrn这里n是能带的指标，k是Bloch波矢，它们),( rkun是r的周期函数，即都是好量子数，可作为波函数的标号，而),(),(rkuRrkunln从另一方面来说，波函数 及其对应的能量本征值 都可看做Bloch波矢的函数，而且满足

33、如下的周期性),(rkn)(kEn),()(),(),(kEKkErkrKknmnnmn这里 是波矢空间的倒格矢。因此，波函数和能量都是波矢空间的周期函数，我们只需计算出一个周期区域中的波函数和能量就够了，这个周期区域称为布里渊区。mKl , m , n 分别是两原子连线在x, y, z上的方向余弦。原子轨道间的哈密顿矩阵元 （一） 离子赝势 对原子的全电子势问题，总的电子数密度 可以分成芯电子数密度 和价电子数密度 之和，即( ) r( )cr( )Vr（二） 模型赝势 可以看出，离子赝势的形式是非常复杂的，而最简单的离子赝势形式是定域的、且与能量无关的：2( )( )()( )( )2(

34、,)( ,)( )PSPSPSVijjijiijjPScCoulcxccjZVrrTrrrZVrVrrr 这里 作为一个可调参数用来拟合原子数据。另一相似的形式是cr这种离子赝势用动量空间变量 表示的形式近似为q据此，加入一些可调参数可以构造模型赝势，称为半经验赝势。它的一般形式为加入 因子使我们可以通过适当地选择( )以保证这种赝势的傅里叶展开很快收敛。常用的一些半导体材料的离子膀势的参数久可从文献上得到。这种形式的赝势通常称为软芯(softcore)势。上述形式的赝势也是定域的。一般地离子赝势应是非定域的，但上述形式已是一个可行的较好的近似。44exp()a q40a 芯态电子价电子赝势和

35、赝波函数只能描写价电子在核间区的行为赝势和赝波函数只能描写价电子在核间区的行为,而不能描写芯态电子的行为。而不能描写芯态电子的行为。 赝势的自洽计算不能包括芯态电子的密度。赝势的自洽计算不能包括芯态电子的密度。 芯态电子密度的自洽只能在孤立原子内进行。芯态电子密度的自洽只能在孤立原子内进行。？如何在全固体的赝势计算和孤立原子的芯态如何在全固体的赝势计算和孤立原子的芯态计算之间进行区分和连接？计算之间进行区分和连接？ 如果所有的电子全被处理成芯电子，就没有屏蔽，则2( )( )ionllVrVrrrr 处的电场：2204( ) ( )rE rdr rrrr 处的库仑势：( )( )rCoulom

36、bVrdr E r 分部积分：2202020201( )4 ( )14 ( )4 ( )4 ( )4 ( )4 ( )rrCoulombrrrrrrVrdrdr rrrddr rrrdr rrdr rrrdr rrdr rrr 下图为产生从头算离子赝势方法的示意图。在芯区外，赝势、赝波函数与真实势、真实波函数完全复合。 孤立原子孤立原子 全电子势全电子势自洽计算孤立原子全电自洽计算孤立原子全电子波函数和全电子势子波函数和全电子势 可调参数：芯半径可调参数：芯半径构造截断函数构造截断函数 抹平势阱抹平势阱 使本征值一致，且使本征值一致，且波函数在芯半径外波函数在芯半径外除比例因子外一致除比例因子

37、外一致构造过渡原子赝势构造过渡原子赝势 归一化条件归一化条件 截断区加上截断区加上修正项修正项模守恒的赝波模守恒的赝波函数函数 由模守恒的赝波函数构造构造模守恒赝势构造模守恒赝势 得到孤立离得到孤立离子的赝势子的赝势扣除价电子对扣除价电子对赝势的贡献赝势的贡献 由上述产生从头算赝势过程可以看出，本质上从头算原于赝势是核与芯电子联合产生的有效势，是从原子的薛定锷方程从头计算得到的，这种赝势可以给出价电子或类价电子(包括部分芯电子，如果需要的话)的正确的电荷分布，因此适合作自治计算。它具有较好的传递性，可用在不同的化学环境中但它的定域性较强，使得动量空间的展开收敛较慢。第二步自洽计算： 在整个固体

38、中，整个赝势包括：（1）所有孤立离子的赝势 + （2）价电子所产生的库仑势和交换关 联势。 在自洽计算中，（1）保持不变， 对（2）进行自洽计算。Schrodinger或 KS方程求出能量本征值和本征波函数确定电子的占据情况从本征波函数计算电子密度从电子密度计算电子所产生的势能初始势能比较代入和得到的势能一致不一致输出能量本征值、波函数、势能 （四）模守恒赝势的解析形式 为了应用方便，从头计算得到的模守恒赝势的数值结果还需拟合成解析函数的形式，分定域部分和非定域部分第六节 缀加平面波方法一、势场和基函数的选取 （一）Muffintin势 固体的能带计算实际包括了两部分主要内容：一是建立一个合理

39、的单电子哈密顿量；二是求解薛定锷方程或KohnSham方程包括将晶体波函数用合理的基函数展开。除了前面二章介绍的几种能带计算方法之外，还有一类方法，它们以一个原胞中电子的能量和波函数为出发点，将晶体的波函数用原胞中电子波函数为基函数展开，并建立晶体波函数在原胞边界面上所必须满足的条件，由此来确定晶体波函数或是基函 数中的某些展开系数。从这一思想出发，发展了原胞法、缀加平面波方法(APW)及格林 函数方法等能带计算方法。原胞法目前已很少被应用、但它的一部分思想却被APW等方法加以借鉴。 Slater提出了Muffin-tin势(蛋糕模子)。他的主要思想是把原胞分为两个区域，以原子为中心的球形区及

40、球外的区域。取一个以原子为中心、半径为 的球。关于 的选取，后面再讨论只是要求各个原子球互不相交。在球内，取球对称势；球外则取常数势。通常选取适当的能量零点，使此常数为零。这样的势场模型称为Muffintin势，因为它很象蛋糕模子而得名。在一个原胞中势场可以表示为 Muffin-tin势的选取可以有不同的方法，常用的是如下的取法，对一个原子周围的势场贡献最大的是中心原子的势场，然后还有它的最近邻原子对这部分空间的势场。次近邻和远邻原子的贡献逐渐减弱。 LFMattheiss 提出了一个构成Muffin-tin势的方法，在许多能带计算中取得了成功的结果，尤其对金属体系更好。他的主要思想是在中心原

41、子的势场上叠加上周围原子势场的以中心为原点的球谐函数展开的首项Y00，即相对于原点的球对称部分。 设空间一点P相对于原于1 和原子2 的座标分别为r1和 r2，原子1和2间的距离为a(见图)。将中心位于原子1的函数 展开为以原于2为中心的各级球谐函数：只保留了lm0的项，即中心原子库仑项的贡献为根据原子自治场计算求得，然后再解泊松方程求U0(r)。相邻其他原子的贡献则根据上面所说的方法叠加上去，总的库仑势为其中 (r)为晶体的空间电荷密度，计算时的切始数值是由中心原子的电荷密度加上相邻原子的电荷密度来求的：（二）缀加平面波 将原胞分为两个区：球内区域称为区域I，球外为区域II。在区域I 有球对

42、称势，KohnSham方程的解应有如下形式：在第 个球内，APW函数可以定义为各个分波函数的线性组合：在球外，势场为零、解应有平面波的形式。设第 个球球心的位矢为 ，则 ，所以在第 个球外，rrr 在球表面上，根据波函数连续的条件，可以得到关于系数 的方程：lma于是，APW基函数可以写为 这里并没有要求波函数导数在球面上连续，不过这不连续性只是在单个基函数中存在，由这些基函数线性组合构成的晶体波函数仍可以是一个连续并有连续导数的函数。（三）久期方程哈密顿：建立能量变分方程： 能量：3*3*()()IIIIIIIIIIIIsphereIIIIIIEd rd rHdS 已考虑球面的导数不连续。令

43、保留一级小量，tItIIItIIEEE3*3*3*3*3*()()ttttItIIIIIIIItItIIIIItIIIIIItsphereEd rEd rd rd rHd rHdS 由变分原理，有考虑基函数在球面连续，得到能量变分方程用APW基函数展开式代入，得令3*3*1()2IIIIIIIIIIIISEd rd rHdS *3*3*,*,1()2MMijijijiji ji jIIIIIIMIIIIIIijiijji jSEc cd rc cd rHc cdS 得由变分确定系数，二、APW方程 （一）APW矩阵元 前面得到的久期方程，应先算出各个矩阵元。为方便计，将与体积分有关的矩阵元分为一区和二区两部分： ,III