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

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