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

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