《DFT-10复杂体系的O(N)算法【优选资料】》由会员分享，可在线阅读，更多相关《DFT-10复杂体系的O(N)算法【优选资料】（27页珍藏版）》请在金锄头文库上搜索。
1、第第10章章 复杂体系的复杂体系的O(N)算法算法1。引言。引言2。O(N)算法的物理基础算法的物理基础 量子力学局域性量子力学局域性3。O(N)算法的基本策略算法的基本策略4。DFT框架下的框架下的O(N)算法算法5。计算流程和主要步骤。计算流程和主要步骤1应用21。引言。引言Order-N算法或算法或O(N)算法的必要性算法的必要性 目前，目前，DFT第一性原理计算方法，如第一性原理计算方法，如fplapw, fplmto, Car-Parrinello, 从头赝势以及许多量子化学计算方法，对于由大从头赝势以及许多量子化学计算方法，对于由大量原子组成的复杂体系已经不能满足需要。量原子组成的
2、复杂体系已经不能满足需要。 原因是以上传统方法的原因是以上传统方法的数值运算工作量（操作数）数值运算工作量（操作数）Nat3。即体系的原子数增加一倍，必须消耗即体系的原子数增加一倍，必须消耗8倍倍cpu时间。时间。研究计算操作数与体系原子数成比例的方法，即研究计算操作数与体系原子数成比例的方法，即O(N)算法对算法对于研究复杂体系十分迫切。于研究复杂体系十分迫切。本章着重与分析本章着重与分析O(N)算法的算法的物理基础物理基础、实现、实现O(N)算法的算法的基基本策略本策略以及把以及把O(N)算法纳入算法纳入DFT框架框架的方法。的方法。2应用22。O(N)算法的物理基础算法的物理基础量子力学
3、局域性量子力学局域性1.Kohn的的“近视原理近视原理(near-sightedness principle)” （Kohn, PRL 76, 3168 （1996） Kohn证明了如下原理：证明了如下原理：多电子体系的某部分的物理性质，多电子体系的某部分的物理性质，不因远离它的区域有势的变化而受影响不因远离它的区域有势的变化而受影响。 rv(r)r考虑一个量子多粒子系，在考虑一个量子多粒子系，在r处的静态物理性质为处的静态物理性质为F，它依赖，它依赖于于r周围线度为周围线度为 的体积内的坐标，的体积内的坐标， 为为de Broglie波长量级。波长量级。Kohn证明，证明，F对于对于r处势的
4、变化处势的变化v(r)是不敏感的。所以，是不敏感的。所以，r处的势保持不变，但比处的势保持不变，但比 更远的区域会变。更远的区域会变。v(r)=03应用2量子力学局域性量子力学局域性例：大多数量子力学静态性质有局域性：例：大多数量子力学静态性质有局域性：1.分子或固体中的化学键分子或固体中的化学键2.局域态密度局域态密度3.电荷密度分布电荷密度分布4.局域磁矩局域磁矩5.结合能，。结合能，。 它们都只依赖于几个近邻原子它们都只依赖于几个近邻原子“壳层壳层”内的局域环境。内的局域环境。2。局域性的描述。局域性的描述 主要方案：采用局域化的主要方案：采用局域化的Wannier函数和密度矩阵方法。函
5、数和密度矩阵方法。 Wannier函数的衰减行为：函数的衰减行为： 有带隙的绝缘体（有带隙的绝缘体（1D, 3D, 无序，团簇，缺陷和表面），都有无序，团簇，缺陷和表面），都有 指数衰减行为指数衰减行为4应用2局域性的描述局域性的描述一般采用广义一般采用广义Wannier函数（函数（GWF,wi, 非周期系非周期系的局域的局域Wannier函数）构造密度矩阵：函数）构造密度矩阵：N是体系每个自旋的电子数。因为是体系每个自旋的电子数。因为wi是局域化的，是局域化的， 将按将按|r-r|衰衰减。对于绝缘体和金属，减。对于绝缘体和金属， |r-r|都表现出指数衰减率。都表现出指数衰减率。T 0时，时
6、，衰减甚至更迅速。衰减甚至更迅速。在在DFT下，核心问题是使下，核心问题是使 成为一个投影算符，其作用是把它成为一个投影算符，其作用是把它投射到占有态空间。投射到占有态空间。这在数学上等价于要求这在数学上等价于要求 必须是等幂的必须是等幂的 （Idempotent）,即要求其本征值在（即要求其本征值在（0,1）区间。）区间。如何把一个接近等幂的密度矩阵如何把一个接近等幂的密度矩阵 变为等幂矩阵变为等幂矩阵 将在下面将在下面介绍。介绍。(10.1)5应用23。O(N)算法的基本策略算法的基本策略如何实现如何实现O(N)算法？算法？vivi根据根据Kohn近视原理，把体积为近视原理，把体积为V的体
7、系分成的体系分成N个子体积个子体积vi(i-1.N)vi 3. 在在vi处取体积为处取体积为vi的区域，它包括的区域，它包括vi和一个缓冲区，然后和一个缓冲区，然后解出每一个解出每一个vi的性质。如果的性质。如果vivi，那么，那么vi内的性质是相当精确内的性质是相当精确的。由于计算每一个的。由于计算每一个vi的工作量完全独立于体系的大小，只要知的工作量完全独立于体系的大小，只要知道道vi内的资料即可。整个体系的大小内的资料即可。整个体系的大小vi的数目的数目N, 于是得到线性于是得到线性标度算法。标度算法。V=N*vi6应用2O(N)算法的基本策略算法的基本策略考虑到处理波函数和密度矩阵的具
8、体要求，考虑到处理波函数和密度矩阵的具体要求，已经提出了多种已经提出了多种O(N)算法方案：算法方案：1.Fermi算符展开方法（算符展开方法（FOE）2.Fermi算符投影方法（算符投影方法（FOP）3.分治（分治（Divide and conquer, D&C）方法）方法4.密度矩阵最小化方法（密度矩阵最小化方法（DMM）5.轨道最小化方法（轨道最小化方法（OM）6.优化基密度矩阵最小化方法（优化基密度矩阵最小化方法（OBDMM）采用采用Chebyshev多多项式将项式将DM展开展开杨伟涛教授杨伟涛教授与与DFT密切结合密切结合7应用2McWeeny净化算法净化算法McWeeny提出一种将
9、接近等幂的密度矩阵提出一种将接近等幂的密度矩阵 变换为更接变换为更接近等幂的密度矩阵近等幂的密度矩阵 的算法：的算法：用用 和和 分别表示分别表示 和和 的本征值，这两个本征值的关系的本征值，这两个本征值的关系是是可见可见所以，这种映射迭代将驱使本征值趋于所以，这种映射迭代将驱使本征值趋于0或或1，由此得到符合，由此得到符合等幂要求的等幂要求的 。(10.2)(10.3)8应用2LNV密度矩阵最小化方法密度矩阵最小化方法Li, Nunes, Vanderbilt (LNV)提出提出DMM方法，文方法，文献上常称献上常称LNV方法。它所采用的净化方法有完全方法。它所采用的净化方法有完全不同的方式
10、。其线性标度是通过对密度矩阵的空不同的方式。其线性标度是通过对密度矩阵的空间截断得到的。间截断得到的。 Ref. Li, Nunes, Vanderbilt : Phys. Rev. B47, 10891 (1993)LNV方法已经在方法已经在TB方法的框架下得到广泛应用。方法的框架下得到广泛应用。采用化学势固定使总能最小的方法。采用化学势固定使总能最小的方法。(后来发现，后来发现，固定化学势方法在实际计算上并不是最方便的固定化学势方法在实际计算上并不是最方便的)。 Ref. Nunes, Vanderbilt : Phys. Rev. B50, 17611 (1994)9应用2线性标度的线性
11、标度的HGG方法方法-1HGG（Hernndez-Gillan-Goringe）方法属于自）方法属于自洽第一性原理方法，并与洽第一性原理方法，并与LNV方法密切相关。方法密切相关。方法的特点方法的特点：1.基态的描述：把基态的描述：把DFT中关于总能中关于总能Etot是是KS占有轨占有轨道道 i 或电子密度的泛函，等价地表述为密度矩阵或电子密度的泛函，等价地表述为密度矩阵的泛函。并要求密度矩阵是等幂的。的泛函。并要求密度矩阵是等幂的。2.采用局域化采用局域化支持函数（支持函数（support function） i 和空和空间有限的有限的变分参数矩分参数矩阵Li j 来表示密度矩来表示密度矩阵
12、。3.用用变分法求分法求总能关于支持函数和能关于支持函数和Li j 均均为最小。最小。HGG方法采用的是固定方法采用的是固定电子数而不是固定化学子数而不是固定化学势。计算上更算上更为方便。方便。10应用2线性标度的线性标度的HGG方法方法-2用用KS占有轨道定义密度矩阵占有轨道定义密度矩阵由由Etot关于关于 最小化求基态，条件是最小化求基态，条件是 (r,r)为等幂及电子数固为等幂及电子数固定。即定。即可以应用可以应用McWeeny净化方法，使密度矩阵达到等幂要求。净化方法，使密度矩阵达到等幂要求。对于实际的第一原理计算，初始的对于实际的第一原理计算，初始的 必须做成可分离必须做成可分离形式
13、，形式，HGG用支持函数用支持函数 i 和局域变分参数和局域变分参数Li j 表示为表示为和和(10.4)(10.5)(10.6)(10.7)11应用2线性标度的线性标度的HG方法方法-3净化之后称为等幂的密度矩阵净化之后称为等幂的密度矩阵上式矩阵上式矩阵K与与L的关系是：的关系是： K = 3LSL - 2LSLSLS 是交叠矩阵是交叠矩阵为了实现线性标度算法，要求为了实现线性标度算法，要求:1。支持函数。支持函数 i 0, 只在某局域空间范围只在某局域空间范围(称为支持区称为支持区)之内之内。2。Li j 0，只有当相应的区域以，只有当相应的区域以 截断距离截断距离Rcut被分离时。被分离
14、时。由于密度矩阵的衰减行为上述条件一般都能满足。由于密度矩阵的衰减行为上述条件一般都能满足。(10.8)(10.9)(10.10)12应用2线性标度的线性标度的HG方法方法-4以下的步骤就是采用变分法，在电子数固以下的步骤就是采用变分法，在电子数固定的条件下，求总能关于支持函数和定的条件下，求总能关于支持函数和L矩阵矩阵为最小，由此得到真正基态能量的上限。为最小，由此得到真正基态能量的上限。目前的目前的O(N)算法，仅限于基态性质的研究。算法，仅限于基态性质的研究。13应用24。DFT框架下的框架下的O(N)算法算法以上基本原理的实际执行，可以在以上基本原理的实际执行，可以在LDA近似下采用近
15、似下采用赝势法。但是，是在实空间的网格点上计算。赝势法。但是，是在实空间的网格点上计算。以每一个原子为中心，取半径为以每一个原子为中心，取半径为Rreg的球作为支持区。的球作为支持区。每一个支持区包含一定数量每一个支持区包含一定数量v的支持函数，并且各区的支持函数，并且各区的的v都一样。都一样。实际执行表明，支持函数的总数实际执行表明，支持函数的总数 0.5Nel（val）。在原来的方法中，每一个支持函数在原来的方法中，每一个支持函数 i (r)都用它在该都用它在该区的网格点区的网格点rl 上之值上之值 i (rl)表示。后来发现这一方法表示。后来发现这一方法在动能计算精度及不同的网格点总能出
16、现不连续性在动能计算精度及不同的网格点总能出现不连续性等问题。新方法中采用一种局域函数将等问题。新方法中采用一种局域函数将 i 展开。展开。14应用2支持函数的表达式支持函数的表达式支持函数用局域化的基函数展开，他们称这种基函数为支持函数用局域化的基函数展开，他们称这种基函数为“视点函数（视点函数（blip function)”。 Rin是第是第i个原子的支持区内的视点网格点（个原子的支持区内的视点网格点（blip-grid）的位置。）的位置。 在实际计算中，对在实际计算中，对 i 的变分采用对的变分采用对bi n的变分。的变分。 计算方法中的一个关键部分是对在积分网格上一组计算方法中的一个关
17、键部分是对在积分网格上一组rl点的点的 i (r) 计算。这些计算结果将用于矩阵元的计算。计算。这些计算结果将用于矩阵元的计算。 从从blip-grid上上bi n之值变换为积分网格上之值变换为积分网格上 i (r)之值的效率是之值的效率是 借助于将视点函数写成如下乘积实现的：借助于将视点函数写成如下乘积实现的：其中其中x, y, z 是是r的直角坐标，的直角坐标， (x)被选择用被选择用B-spline工作。工作。(10.11)(10.12)15应用2DFT框架下的框架下的O(N)算法算法-2主要计算公式：主要计算公式：其中，动能（对所有网格点求和）：其中，动能（对所有网格点求和）：其它三个
18、能量全部依赖于其它三个能量全部依赖于rl 点上的电子密度点上的电子密度具体计算并不涉及特殊技巧，例如具体计算并不涉及特殊技巧，例如LDA交换关联能可对交换关联能可对求和得到。求和得到。(10.13)(10.14)(10.15)16应用2DFT框架下的框架下的O(N)算法算法-3通过变分法使总能最小，在通过变分法使总能最小，在HGG方法中，采用共方法中，采用共轭梯度近似，涉及如下解析式：轭梯度近似，涉及如下解析式：(10.16)(10.17)(10.18)(10.19)17应用2以及(10.20)(10.21)(10.22)(1023)18应用2矩阵乘积中矩阵乘积中Hi j 是支持函数是支持函数
19、 i 和和 j 之间之间的的KS-Hamiltonian矩阵元。矩阵元。线性标度来自支持函数的空间局域性。因线性标度来自支持函数的空间局域性。因为支持函数之间的距离为支持函数之间的距离超过某一截断距离时，超过某一截断距离时，H和和S都将都将0。19应用2DFT框架下的框架下的O(N)算法算法-4用以上截断值到矩阵用以上截断值到矩阵L上，上，Etot及其微商的及其微商的表达式中的所有矩阵都是稀疏矩阵。表达式中的所有矩阵都是稀疏矩阵。稀疏矩阵的非稀疏矩阵的非0矩阵元的数目矩阵元的数目 原子数（成原子数（成线性比例）。由此得到线性标度的算法。线性比例）。由此得到线性标度的算法。O(N)算法正成为研究
20、大原子数复杂体系的算法正成为研究大原子数复杂体系的有力工具。有力工具。20应用25。计算流程和主要步骤。计算流程和主要步骤1.计算流程图计算流程图2.主要计算步骤主要计算步骤21应用21. General computational strategy.计算计算E和和dE/dL选择选择L空间的搜索方向空间的搜索方向E关于关于L最小化计算最小化计算检查检查E关于关于L是否收敛是否收敛检查检查E关于关于 是否收敛是否收敛End计算计算dE/d , 修改修改 输入猜测的输入猜测的L和和 内循环内循环外循环外循环NoNo22应用22. 主要计算步骤主要计算步骤1.在在O(N)算法中有两个变量：算法中有两
21、个变量： 和和L。2.Etot对对L的变分是在的变分是在 固定条件下进行的，它固定条件下进行的，它构成计算的构成计算的内循环内循环。要求。要求Etot最小，这需最小，这需要在要在L空间搜索一系列方向，每一个搜索的空间搜索一系列方向，每一个搜索的方向有方向有dE/dL决定。决定。3.Etot对对 的变分是的变分是外循环外循环。改变。改变 使使Etot最小。最小。23应用2内循环的步骤内循环的步骤每一次改变每一次改变L由如下步骤组成：由如下步骤组成：1.由方程由方程(10.16)-(10.19)计算计算 Etot/ L，并用，并用它计算新的它计算新的L。2.由新的由新的L计算计算K（10.9）。）
22、。3.利用利用K计算计算Etot。4.如果如果Etot的变化小于设定的的变化小于设定的tolerance，循环，循环终止，否则，利用终止，否则，利用n(r)重新计算重新计算H矩阵（矩阵（S矩阵不变，因为在内循环矩阵不变，因为在内循环 i 是固定的），是固定的），然后回到第一步。然后回到第一步。24应用2外循环的步骤外循环的步骤1.由方程由方程(10.20)-(10.23)计算计算 Etot/ bi n。并用它计算新的并用它计算新的bi n。2.用新的用新的 计算计算S，并计算新的，并计算新的K。3.用用K计算新的计算新的Etot。4.判断是否满足收敛条件，不满足时，修改判断是否满足收敛条件，不满足时，修改 ，进行第，进行第2步。步。5.计算新的计算新的KS势和势和H。6.如果如果Etot的变化小于设定的的变化小于设定的tolerance，循，循环终止，否则，回到第环终止，否则，回到第1步。步。25应用2阅读参考阅读参考26应用2End27应用2
