计算物理方法(sec1, sec2)

计算物理熊诗杰（用于凝聚态物理）§ 凝聚态物理中数值计算的基本特点 § 大型矩阵的数值严格对角化方法 § 密度泛函理论和固体能带理论 § Monte Carlo方法第一章 凝聚态物理中的数值计算一、凝聚态系统基本特点 1. 由大量粒子组成，是一个多粒子体系；2. 这些粒子都是微观粒子，满足量子力学的运动 规律。因此，凝聚态系统的状态要用多体波 函数表达。二、多体波函数满足的基本方程- 薛定谔方程这里， 是多体波函数，其自变量是时间和描写各 个粒子状态的广义坐标， 是系统的哈密顿算符。如果 不依赖时间，有我们有如下的定态薛定谔方程：是能量的本征值。三、薛定谔方程的基本特点1. 是关于未知波函数的线性齐次方程-如果选取一组正交归一完备的基 函数，可以将波函数表达为向量，哈密顿 算符表达为矩阵。2. 只有 取一定的值才有非零的波函 数- 能量本征值3. 是厄米算符，能量本征 值为实数。4. 如果哈密顿表达为 阶矩阵，则有 个能量本征值 ，对应着 个正交归一的本 征波函数。本征值有可能在数值 上一样（简并）。求能量的本征值和对应的本征波函数。波函数可以是单体的波函数，也可以是多体的波 函数。数值计算的基本出发点，是构造一组完备（或近 于完备）的基函数，使波函数成为向量，哈密顿 成为矩阵。四、数值计算的基本任务第二章：大型Hamiltonian系统的 数值严格对角化方法Lanczos 方法第一节 概述 在对一个实际的物理系统进行理论研究的过程中 ，经常会遇到要求一个大型或超大型Hermit矩阵的 本征值和本征函数的问题。矩阵的维数往往会达到 几十万甚至上百万。这样的矩阵在包含有大量微观 自由度的宏观体系中是经常遇到的。例如，如果我们用分立的自旋来描写磁性 晶格中的磁性原子，而磁性原子间的互作 用则用相邻的自旋格点间的交换积分来表 达，我们可以得到如下的Heisenberg哈密 顿量： 这里Jx ,Jy , Jz是三个方向上的耦合常数，Sxi , Syi , Szi 分别是第i个格点上的自旋在x, y, z三个方向 上的分量。如果我们要将这个哈密顿量用数值的 方法对角化，以求出其本征值和本征函数，我们 就必须首先选择Hilbert空间中一组完全的向量集作为基，把哈密顿量在这个基上用矩阵 的形式表达出来，再对该矩阵进行数值求 解。然而，与能带论中所用的单电子近似 不同的是，这种强关联的多粒子体系所对 应的Hilbert空间的基矢的个数是极其庞大 的。在我们这个例子中，每个格点上的自 旋的状态数是(2S+1)，如果我们研究的有 限的系统包含有N个格点，那么总的状态数 就是(2S+1)N个，也就是说，哈密顿矩阵的 大小应该是(2S+1)N×(2S+1)N。这是一个极 其庞大的数字，例如，我们研究一个只 有4×4大小的自旋为1/2的二维晶格，那么 哈密顿矩阵的大小就是65536×65536，这 样大小的矩阵的直接数值求解在现在的绝 大多数计算机中都是办不到的。矩阵维数的增大使得求解本征值和本 征函数变得非常困难，但是，我们看到， 在绝大多数在物理上感到兴趣的问题中， 哈密顿矩阵都有如下特点：(1) 哈密顿矩阵都是厄米矩阵。(2) 尽管哈密顿矩阵异常庞大，但是其中绝大 多数非对角矩阵元都是零，只有数量很少的非 零矩阵元，也就是说，这种矩阵是极为稀疏的 矩阵。例如，对我们刚才讲的例子，哈密顿矩 阵的每一行尽管有65536个元素，但其中非零的 矩阵元不会超过5个。(3) 所考虑的物理系统，一般都具有一定的对 称性，满足一定的守恒律。这样，我们就可以 根据对称性和守恒律的要求，采用一定的基函 数，使哈密顿矩阵分解成为互不关联的子矩阵 ，从而可以对比较小的子矩阵分别求解。Lanczos方法就是根据上述特点所发 展起来的一种极为有效的哈密顿矩阵求解 本征值和本征函数的数值方法。近年来， 它已被广泛应用到强关联系统、无序系统 等极为重要的物理系统之中，成为计算物 理中一项至关重要的新方法。我们这一单 元，就是介绍这一方法的基本做法和它的 一些重要的应用。第二节 基本方法发展Lanczos方法的最重要的目的，就是要找 出一条有效的途径，来解决哈密顿矩阵过于庞大的 问题。问题的关键就在于，如何充分利用前面说到 的三个特点，即厄米性，稀疏性和对称性，来达到 我们的目的。这是一个极其复杂的问题，涉及到一 个物理系统的一些非常具体的特点。在这一节中， 我们将着重介绍如何利用体系的厄米性来将一个庞 大的哈密顿矩阵约化为较小的矩阵，这是对各种物 理系统都通用的基本的Lanczos方法。至于利用稀 疏性和对称性来作进一步的约化，我们将在以下几 节中加以介绍。 我们设H是某一在N维Hilbert空间中的 厄米算符，而v1是该空间中的一个任意的归 一化的向量，即：v1+v1=1。我们连续地用 H去作用，来构造出一组正交归一的向量：.我们注意到这些方程的结构是三对角型的，即 每个方程的右边只出现相邻的三个元素。这是由哈 密顿矩阵的厄米性所决定的。这一迭代序列将在正 交归一向量走遍整个Hilbert空间后自动中止，也就 是说，在做到第N步时，我们有： 但由于在Hilbert空间中最多只能有N个相互正交归一 的向量，因此新的向量vN+1必须为零。这样，这一方 法具有一种十分有用的性质，即：即使我们一开始并 不知道向量空间的维数，这一方法也会很好地自动中 止。注意到我们这里处理的是一般的算符，而不是矩 阵，因此空间的维数不需要明显地写出来。 这些向量，v1, ,vN , 称为Lanczos向量 ，它们构成了Hilbert空间的一组正交归一的 完备集，以它作为基，哈密顿算符可表达 为三对角矩阵的形式：如果U是H的本征向量，而u是 的对应的本征向量 ，它们之间的关系是：UVu这里V是以Lanczos向量作为列的矩阵。如果算符H具有简并的本征值，则初始向量v1将 包含对应于这个简并的本征值的一个本征向量（这 是因为v1在这个简并本征值所对应的子空间的投影 只有一个线性无关的向量）。这意味着Lanczos过 程将在小于N步的地方中止，而且这个过程只给出 对应于这个简并集合中的一个本征值和一个本征向 量。一般来讲，如果H具有n个不同的本征值 （n<N）E1，E2，, En，它们分别是p1, p2, ,pn 重简并，那么上述的迭代过程将中止在第n步，所 得到的由Lanczos向量所张开的n维子空间将只包 含每个简并集合中的一个本征态。所得到的n×n 阶三对角矩阵将具有n个不同的本征值。于是，我 们需要再去求出剩下的N-n个本征值和本征向 量。我们可以选择一个新的初始向量，使它与所 得到的所有Lanczos向量正交，来继续上述过 程。同样，这个过程将经过m步之后中止，m应 等于所剩下的具有不同数值的本征值的个数。我 们可以这样继续下去，直到得到整个Hilbert空间 中的全部本征值和本征向量为止。这样得到的三 对角矩阵具有块对角形的结构，即：这里每个块都是非简并的，而且其维数等于在前 面的块构成之后所剩下的不相等的本征值的个 数。 Lanczos方法具有一个十分有趣而且至关重要 的性质，即：在经过p - 1步后，所得到的p个 Lanczos向量可以完全决定三对角矩阵的左上方 的大小为p×p的部分的元素，如果我们对这个p×p 的矩阵求本征值和本征向量，那么在p增大的时候 ，这些本征值和本征向量将会极快地收敛于完全 的矩阵的低端的本征值和本征向量。对一个物理 系统，在多数情况下，我们只对能量较低的本征 值和本征向量感到兴趣，因为它们对应于系统的 基态和能量较低的激发态。这样，我们就可以根 据Lanczos方法的这个性质，不去管原先的矩阵H 有多大，只要做有限的p步，得到p×p大小的矩阵 ，求出它的本征值和本征向量，然后增大p的值， 使得到的低端的本征值收敛。只要它们收敛，我们就可以认定所得到的低端本征值和本征向量就 是完全矩阵的低端本征值和本征向量，从而得到所求系统 的基态和较低的激发态。上述的Lanczos方法的这一性质，已经得到了严格的 证明。在计算实践中，收敛速度之快令人吃惊。在许多情 况下，不管的大小有多大。所取的Lanczos步数往往只要 取50到100之间，收敛的精度就可以达到计算机的精度。 但这并不是说Lanczos方法很容易做。实际上，Lanczos 方法是计算物理中技巧性非常高的工作。这一方面是因为 所研究的物理对象极为复杂，有时光是为了写出矩阵元就 要费很大的工夫。另一方面，在算法中某些不引人注目的 误差，如舍入误差等，都会破坏整个程序的运行。因此， 上述的Lanczos方法只有在运用精确的算法和程序时才是 有意义的。在实际的计算中，计算误差将以各种方式影响到整个过程。这里我们暂不讨论Lanczos方法 本身的精确性和稳定性，但我们指出一些不精确 的算法所产生的后果，初看起来，这些后果似乎 是十分奇怪的。例如，由于计算过程中产生的舍 入误差，使得在N1步后，所得到的vN+1并不准 确为零，这样整个过程就不会在N1步后中止。 尽管这里的vN+1并不大，但已足以妨碍整个过程 的正常运行。如果继续进行下去，同样的本征值 将会被一次又一次地产生。类似地，虽然从理论 上来说在第N1步中vN必须与前面所得到的所有 向量正交，然而在实践上vN是通过以下的公式计 算的：在本式中由于可能出现两个相近的数字相 减，会引起较大的有效数字的丧失，最后 的结果可能是非常奇异的。标准的方法是 ：为了保证其正交性必须将当前的向量v 与前面的所有的向量从新进行正交化。即 ，必须在程序中引入如下的步骤：§ 这是重正交化手续，它能保证Lanczos方法能正 常中止。 § Lanczos方法使得我们能够把一个很大的矩阵 转化为一个等效的、维数较小的三对角矩阵。为 了求出对应的本征值和本征态，我们需要一个有效的数值方法，以便对三对角矩阵进行对角 化。在下一节，我们将给出两种较为常用的方法 ，对三对角矩阵进行对角化。第三节 三对角矩阵的对角化在这一节中，我们研究两种对三对角矩阵 进行对角化的方法。下面我们将分别对它们进 行阐述。(1) Dean方法。在Lanczos方法中，我们得到了有限维数 的三对角矩阵，由于初始的Hamilton矩阵是 Hermit矩阵，所以得到的三对角矩阵也是 Hermit矩阵。这种三对角的Hermit矩阵可以很 方便地用Dean方法来进行对角化。这种方法 是基于如下的负本征值定理：定理如果矩阵是M阶的三对角Hermit矩阵，则 它的比某一实数 小的本征值的个数等于所有的 实数U(i) (i=1, 2, 3, ., M) 中的负数的个数，这里 实数U (i) 由如下的递推关系来决定：E(i) 是矩阵 的第i个对角矩阵元，而t(i) 是矩阵 的第i个次对角矩阵元。 根据这个定理，我们可以方便地通过改变 的值 求出我们所需要的能量范围里本征值的分布情况，并 且可以求出各个本征值的准确位置。在求出本征值后 ，我们可以通过如下的方法求出所对应的本征函数： 我们设定对应于本征值 的本征函数在Lanczos 向量所张开的空间中的第n个分量是1，则其它的分 量可以通过如下的递推关系求出：最后再通过归一化，即可得到对应的本征函数。 实际上，一旦我们原来的厄米矩阵已通过Lanczos 方法约化为三对角形，有关确定其本征值的可能的 方式就是直接寻找特征多项式pN()的根。上述方法 就是通过递推过程产生出低次多项式的Sturm序列 ，该序列将本征值局限于一些实轴区间上，然后可 以运用象二分法或Newton法那样的求根方法来精炼 区间。对应的本征向量于是也可以用上述的逆迭代 法求得。 （2） 运用标准的Lapack或Linpack子程序计算 三