DFT-10复杂体系的O(N)算法PPT课件

1、1第第10章章 复杂体系的复杂体系的O(N)算法算法1。引言。引言2。O(N)算法的物理基础算法的物理基础 量子力学局域性量子力学局域性3。O(N)算法的基本策略算法的基本策略4。DFT框架下的框架下的O(N)算法算法5。计算流程和主要步骤。计算流程和主要步骤21。引言。引言 Order-N算法或算法或O(N)算法的必要性算法的必要性 目前，目前，DFT第一性原理计算方法，如第一性原理计算方法，如fplapw, fplmto, Car-Parrinello, 从头赝势以及许多量子化学计算方法，对于由大从头赝势以及许多量子化学计算方法，对于由大量原子组成的复杂体系已经不能满足需要。量原子组成的复

2、杂体系已经不能满足需要。 原因是以上传统方法的原因是以上传统方法的数值运算工作量（操作数）数值运算工作量（操作数）Nat3。即体系的原子数增加一倍，必须消耗即体系的原子数增加一倍，必须消耗8倍倍cpu时间。时间。 研究计算操作数与体系原子数成比例的方法，即研究计算操作数与体系原子数成比例的方法，即O(N)算法算法对于研究复杂体系十分迫切。对于研究复杂体系十分迫切。 本章着重与分析本章着重与分析O(N)算法的算法的物理基础物理基础、实现、实现O(N)算法的算法的基基本策略本策略以及把以及把O(N)算法纳入算法纳入DFT框架框架的方法。的方法。32。O(N)算法的物理基础算法的物理基础量子力学局域

3、性量子力学局域性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)=04量子力学局域性量子力学局域性例：大多数量子力学静态性质有局域性：例：大多数量子力学静态性质有局域性：分子或固体中的化学键分子或固体中的化学键局域态密度局域态密度电荷密度分布电荷密度分布局域磁矩局域磁矩结合能，。结合能，。 它们都只依赖于几个近邻原子它们都只依赖于几个近邻原子“壳层壳层”内的局域环境。内的局域环境。2。局域性的描述。局域性的描述 主要方案：采用局域化的主要方案：采用局域化的Wannier函数和密度矩阵方法。函数和密度矩阵方法。 Wannier

5、函数的衰减行为：函数的衰减行为： 有带隙的绝缘体（有带隙的绝缘体（1D, 3D, 无序，团簇，缺陷和表面），都有无序，团簇，缺陷和表面），都有 指数衰减行为指数衰减行为5局域性的描述局域性的描述 一般采用广义一般采用广义Wannier函数（函数（GWF,wi, 非周期系非周期系的局域的局域Wannier函数）构造密度矩阵：函数）构造密度矩阵：*( , )2( )( )Niiir rw r w rN是体系每个自旋的电子数。因为是体系每个自旋的电子数。因为wi是局域化的，是局域化的， 将按将按|r-r|衰衰减。对于绝缘体和金属，减。对于绝缘体和金属， |r-r|都表现出指数衰减率。都表现出指数衰减

6、率。T 0时，时，衰减甚至更迅速。衰减甚至更迅速。在在DFT下，核心问题是使下，核心问题是使 成为一个投影算符，其作用是把它成为一个投影算符，其作用是把它投射到占有态空间。投射到占有态空间。这在数学上等价于要求这在数学上等价于要求 必须是等幂的必须是等幂的 （Idempotent）,即要求其本征值在（即要求其本征值在（0,1）区间。）区间。如何把一个接近等幂的密度矩阵如何把一个接近等幂的密度矩阵 变为等幂矩阵变为等幂矩阵 将在下面将在下面介绍。介绍。(10.1)63。O(N)算法的基本策略算法的基本策略 如何实现如何实现O(N)算法？算法？vivi根据根据Kohn近视原理，把体积为近视原理，把

7、体积为V的体系分成的体系分成N个子体积个子体积vi(i-1.N)vi 3. 在在vi处取体积为处取体积为vi的区域，它包括的区域，它包括vi和一个缓冲区，然后和一个缓冲区，然后解出每一个解出每一个vi的性质。如果的性质。如果vivi，那么，那么vi内的性质是相当精确内的性质是相当精确的。由于计算每一个的。由于计算每一个vi的工作量完全独立于体系的大小，只要知的工作量完全独立于体系的大小，只要知道道vi内的资料即可。整个体系的大小内的资料即可。整个体系的大小vi的数目的数目N, 于是得到线性于是得到线性标度算法。标度算法。V=N*vi7O(N)算法的基本策略算法的基本策略考虑到处理波函数和密度矩

8、阵的具体要求，考虑到处理波函数和密度矩阵的具体要求，已经提出了多种已经提出了多种O(N)算法方案：算法方案：Fermi算符展开方法（算符展开方法（FOE）Fermi算符投影方法（算符投影方法（FOP）分治（分治（Divide and conquer, D&C）方法）方法密度矩阵最小化方法（密度矩阵最小化方法（DMM）轨道最小化方法（轨道最小化方法（OM）1. 优化基密度矩阵最小化方法（优化基密度矩阵最小化方法（OBDMM）采用采用Chebyshev多多项式将项式将DM展开展开杨伟涛教授杨伟涛教授与与DFT密切结合密切结合8McWeeny净化算法净化算法 McWeeny提出一种将接近等幂的密度矩

9、阵提出一种将接近等幂的密度矩阵 变换为更变换为更接近等幂的密度矩阵接近等幂的密度矩阵 的算法：的算法：2332用用 和和 分别表示分别表示 和和 的本征值，这两个本征值的关系的本征值，这两个本征值的关系是是2332可见可见31221212, 0,1 for 0,1 have for 所以，这种映射迭代将驱使本征值趋于所以，这种映射迭代将驱使本征值趋于0或或1，由此得到符合，由此得到符合等幂要求的等幂要求的 。(10.2)(10.3)9LNV密度矩阵最小化方法密度矩阵最小化方法 Li, Nunes, Vanderbilt (LNV)提出提出DMM方法，文方法，文献上常称献上常称LNV方法。它所采

10、用的净化方法有完全方法。它所采用的净化方法有完全不同的方式。其线性标度是通过对密度矩阵的空不同的方式。其线性标度是通过对密度矩阵的空间截断得到的。间截断得到的。 Ref. Li, Nunes, Vanderbilt : Phys. Rev. B47, 10891 (1993) LNV方法已经在方法已经在TB方法的框架下得到广泛应用。方法的框架下得到广泛应用。采用化学势固定使总能最小的方法。采用化学势固定使总能最小的方法。(后来发现，后来发现，固定化学势方法在实际计算上并不是最方便的固定化学势方法在实际计算上并不是最方便的)。 Ref. Nunes, Vanderbilt : Phys. Rev

11、. B50, 17611 (1994)10线性标度的线性标度的HGG方法方法-1HGG（Hernndez-Gillan-Goringe）方法属于）方法属于自洽第一性原理方法，并与自洽第一性原理方法，并与LNV方法密切相关。方法密切相关。方法的特点方法的特点：基态的描述：把基态的描述：把DFT中关于总能中关于总能Etot是是KS占有轨占有轨道道 i 或电子密度的泛函，等价地表述为密度矩阵或电子密度的泛函，等价地表述为密度矩阵的泛函。并要求密度矩阵是等幂的。的泛函。并要求密度矩阵是等幂的。采用局域化采用局域化支持函数（支持函数（support function） i 和和空间有限的空间有限的变分参

12、数矩阵变分参数矩阵Li j 来表示密度矩阵。来表示密度矩阵。1.用变分法求总能关于支持函数和用变分法求总能关于支持函数和Li j 均为最小。均为最小。HGG方法采用的是固定电子数而不是固定化学势。方法采用的是固定电子数而不是固定化学势。计算上更为方便。计算上更为方便。11线性标度的线性标度的HGG方法方法-2 用用KS占有轨道定义密度矩阵占有轨道定义密度矩阵*( , )( )( )iii occr rrr由由Etot关于关于 最小化求基态，条件是最小化求基态，条件是 (r,r)为等幂及电子数固为等幂及电子数固定。即定。即( , ) ( , ) ( , )2( , )elr rdrr rrrNd

13、rr r可以应用可以应用McWeeny净化方法，使密度矩阵达到等幂要求。净化方法，使密度矩阵达到等幂要求。对于实际的第一原理计算，初始的对于实际的第一原理计算，初始的 必须做成可分离必须做成可分离形式，形式，HGG用支持函数用支持函数 i 和局域变分参数和局域变分参数Li j 表示为表示为( , )r r,( , )( )( )iijjijr rr Lr 和和(10.4)(10.5)(10.6)(10.7)12线性标度的线性标度的HG方法方法-3 净化之后称为等幂的密度矩阵净化之后称为等幂的密度矩阵,( , )( )( )iijjijr rr Kr 上式矩阵上式矩阵K与与L的关系是：的关系是：

14、 K = 3LSL - 2LSLSLS 是交叠矩阵是交叠矩阵( )( )ijijSdrrr 为了实现线性标度算法，要求为了实现线性标度算法，要求:1。支持函数。支持函数 i 0, 只在某局域空间范围只在某局域空间范围(称为支持区称为支持区)之内之内。2。Li j 0，只有当相应的区域以，只有当相应的区域以 截断距离截断距离Rcut被分离时。被分离时。由于密度矩阵的衰减行为上述条件一般都能满足。由于密度矩阵的衰减行为上述条件一般都能满足。(10.8)(10.9)(10.10)13线性标度的线性标度的HG方法方法-4 以下的步骤就是采用变分法，在电子数固以下的步骤就是采用变分法，在电子数固定的条件

15、下，求总能关于支持函数和定的条件下，求总能关于支持函数和L矩阵矩阵为最小，由此得到真正基态能量的上限。为最小，由此得到真正基态能量的上限。 目前的目前的O(N)算法，仅限于基态性质的研究。算法，仅限于基态性质的研究。144。DFT框架下的框架下的O(N)算法算法 以上基本原理的实际执行，可以在以上基本原理的实际执行，可以在LDA近似下采用近似下采用赝势法。但是，是在实空间的网格点上计算。赝势法。但是，是在实空间的网格点上计算。 以每一个原子为中心，取半径为以每一个原子为中心，取半径为Rreg的球作为支持区。的球作为支持区。每一个支持区包含一定数量每一个支持区包含一定数量v的支持函数，并且各区的

16、支持函数，并且各区的的v都一样。都一样。 实际执行表明，支持函数的总数实际执行表明，支持函数的总数 0.5Nel（val）。 在原来的方法中，每一个支持函数在原来的方法中，每一个支持函数 i (r)都用它在该都用它在该区的网格点区的网格点rl 上之值上之值 i (rl)表示。后来发现这一方法表示。后来发现这一方法在动能计算精度及不同的网格点总能出现不连续性在动能计算精度及不同的网格点总能出现不连续性等问题。新方法中采用一种局域函数将等问题。新方法中采用一种局域函数将 i 展开。展开。15支持函数的表达式支持函数的表达式 支持函数用局域化的基函数展开，他们称这种基函数为支持函数用局域化的基函数展

17、开，他们称这种基函数为“视点函数（视点函数（blip function)”。( )()ii ninnrbrR Rin是第是第i个原子的支持区内的视点网格点（个原子的支持区内的视点网格点（blip-grid）的位置。）的位置。 在实际计算中，对在实际计算中，对 i 的变分采用对的变分采用对bi n的变分。的变分。 计算方法中的一个关键部分是对在积分网格上一组计算方法中的一个关键部分是对在积分网格上一组rl点的点的 i (r) 计算。这些计算结果将用于矩阵元的计算。计算。这些计算结果将用于矩阵元的计算。 从从blip-grid上上bi n之值变换为积分网格上之值变换为积分网格上 i (r)之值的效

18、率是之值的效率是 借助于将视点函数写成如下乘积实现的：借助于将视点函数写成如下乘积实现的：( )( ) ( ) ( )rxyz其中其中x, y, z 是是r的直角坐标，的直角坐标， (x)被选择用被选择用B-spline工作。工作。(10.11)(10.12)16DFT框架下的框架下的O(N)算法算法-2 主要计算公式：主要计算公式：totkpsHxcEEEEE其中，动能（对所有网格点求和）：其中，动能（对所有网格点求和）：212,2( )()( )kjijriijEdrr Kr 其它三个能量全部依赖于其它三个能量全部依赖于rl 点上的电子密度点上的电子密度,( )2( )( )lilijjl

19、ijn rr Kr 具体计算并不涉及特殊技巧，例如具体计算并不涉及特殊技巧，例如LDA交换关联能可对交换关联能可对( ) ( )lxcln rn r求和得到。求和得到。(10.13)(10.14)(10.15)17DFT框架下的框架下的O(N)算法算法-3 通过变分法使总能最小，在通过变分法使总能最小，在HGG方法中，采用共方法中，采用共轭梯度近似，涉及如下解析式：轭梯度近似，涉及如下解析式：2123()2()()totijelijEDLNELDSLHHLSSLSLHSLHLSHLSLSESLSSLSLS (10.16)(10.17)(10.18)(10.19)18 以及int4()( )(

20、)3()2()tottottoti mi mi mexactgridtotlimilii mgridilijKSijjjijijijEEEbbbErRQrbQrKHGGLHLLSLHLLHLSL (10.20)(10.21)(10.22)(1023)19 矩阵乘积中矩阵乘积中Hi j 是支持函数是支持函数 i 和和 j 之间的之间的KS-Hamiltonian矩阵元。矩阵元。 线性标度来自支持函数的空间局域性。因线性标度来自支持函数的空间局域性。因为支持函数之间的距离为支持函数之间的距离 超过某一截断距离时，超过某一截断距离时，H和和S都将都将0。20DFT框架下的框架下的O(N)算法算法-4

21、 用以上截断值到矩阵用以上截断值到矩阵L上，上，Etot及其微商的及其微商的表达式中的所有矩阵都是稀疏矩阵。表达式中的所有矩阵都是稀疏矩阵。 稀疏矩阵的非稀疏矩阵的非0矩阵元的数目矩阵元的数目 原子数（成原子数（成线性比例）。由此得到线性标度的算法。线性比例）。由此得到线性标度的算法。 O(N)算法正成为研究大原子数复杂体系的算法正成为研究大原子数复杂体系的有力工具。有力工具。215。计算流程和主要步骤。计算流程和主要步骤计算流程图计算流程图主要计算步骤主要计算步骤221. General computational strategy.计算计算E和和dE/dL选择选择L空间的搜索方向空间的搜索方向E关于关于L最小化计算最小化计算检查检查E关于关于L是否收敛是否收敛检查检查E关于关于 是否收敛是否收敛End计算计算dE/d , 修改修改 输入猜测的输入猜测的L和和 内循环内循环外循环外循环NoNo232. 主要计算步骤主要计