电子结构算法.pdf

电子结构计算主要算法 李震宇 (USTC) Outline  单电子方程的离散化 实空间网格  原子基组  平面波基组   赝势 模守恒赝势、超软赝势  投影缀加波方法   自洽场 求解单电子波函数  求解电荷密度   几何优化 http://staff.ustc.edu.cn/~zyli/teaching.html 数值离散  单电子方程 [ 2m 2  vext (r )  vH (r )  v XC (r )]i   ii  实空间网格  边界条件 开放边界条件  周期性边界条件   矩阵形式(本征值问题) 有限差分  中心差分  高阶差分 原子基组  原子轨道的线性组合(LCAO)  (r ,  ,  )  Rn ( r )Ylm ( ,  )  Slater-type orbital (STO)   r, ,   r n1e rYlm 计算三中心和四中心双电子积分比较困难，一般采用数值 积分方法(ADF)  只需要确定Slater指数就可以得到一个基函数，因此对重金 属元素很容易得到可用的基组  Gaussian-type orbital (GTO)  在Gaussian函数前面加上不同的因子  ( x, y, z)  x y z e i j k  r 2 i、j、k之和为0、1、2分别对应s、p、d轨道  5个正则d轨道：xy, yz, xz, x2-y2, 3z2-r2   GTO乘积定理：两个高斯函数的积仍是高斯函数  三维积分变成独立的三个一维积分  GTO并不是一种轨道：r大的时候衰减太快, 而在原子 核处又没有尖点(cusp)。  Gaussian函数没有节点：没法只用一个高斯函数来表 示价轨道 Gaussian基组  收缩高斯基组(STO-3G)  分裂基组(multi-ζ)与分裂价基(3-21G)  极化基组(6-31G**, 6-311G(3df,3pd))  弥散函数(6-31++G)  练习：用6-311++G(3df, 2pd)计算H2O时，用到多少个 轨道，多少个GTO？ 数值原子基组  直接作(赝)原子轨道计算，存储径向函数，样条插值 分子解离成原子exactly  Confining Potential, energy shift, 局域基组   极化、扩展等轨道可以通过离子、激发态、类氢原子、 电场微扰等方法产生。 傅立叶变换  傅立叶变换 F (k )     f ( x )e 2 ikx dx f ( x)     F ( k )e 2 ikx dk  通用的完备基组(平面波基组)  n  r    Cn (g)eigr dg  离散傅立叶变换 N 1 1 F ( k )   f ( x )e N x 0  2  ikx N N 1 f ( x)   F (k )e 2 ikx N k 0 当 f为实数时，有N 个实空间变量，对应倒空间中[1,(N-1)/2] 区间内的,(N-1)/2个复数，以及0和N/2 两个实数。 实空间与k空间对应 实空间与k空间对应 实空间与k空间对应 实空间与k空间对应 实空间与k空间对应 平面波基组  单粒子波函数(非周期)的离散化  周期性边界条件 (布洛赫定理) nk r   unk (r)eikr nk  r    Cn,k G eik G r G  波恩－冯卡曼边界条件 n   N  n  e ikNR 2 1  k  m NR 实空间与倒易空间 哈密顿算符  动能  势能  交换关联势 Evaluate in real space, then FFT to k-space Hartree势 Poisson equation in k-space  Add all contribution, then FFT back to real space  卷积  卷积的定义 f ( x)  g ( x)     f ( ) g ( x   )d 卷积与截断半径  卷积定理 实空间中的乘积对应倒空间中的卷积，两者互为Fourier变 换，反之亦然 1 M 1 f (k )  g (k )  f (m) g (k  m)   倒空间卷积 M m0  Wrap around error  电荷密度(2Gcut)  有效势(2Gcut内有意义) G | veff | G   veff (G n ) GG ,Gn n  单电子方程 赝势  芯电子与价电子 原子形成分子与固体时，芯电子基本保持不变  泡利不相容原理→原子轨道相互正交→价电子态快速振荡   赝势：对赝波函数的散射性质与原子核及芯电子对价 波函数的散射性质相同  赝波函数通常变化缓慢  需要的基函数较少  没有径向节点  无需求解内层电子 赝势构造步骤 原子全电子计算，只需考虑径向部分  (r )  [ul ( r ) / r ]Ylm 2. 构造赝价波函数  ps (r ) 3. 反演得到总的赝势 1.   d 2 ps u r 2  l (l  1) dr 2 l    tot vl ( r )   l   2   2me  r ulps  r     4. 去屏蔽 vlPS  vltot  vH [0PS ]  v xc [ 0PS ]  芯修正 v xc [0PS ]  v xc [ core  0PS ]  部分芯修正 rnlc通常选核密度降到 低于价电子密度的点 模守恒赝势  赝波函数限制条件 光滑，无节点  对某种典型的电子组态，赝本征值与真实本征值相等  真实波函数与赝波函数在芯半径rc外相等  在芯半径内与真实波函数电荷密度积分相同(模守恒)   模守恒条件保证了径向波函数的对数微分对能量的导 数在参考能量处相等  芯半径的选择 Big enough to make soft Psp  Small enough to keep good transferability  Not too small to be very close to the outmost radial node  Kleinman-Bylander近似  不同的角动量感受到不同的势 vˆ PS  | Ylm vl (r )Ylm | lm  在较远处都趋于库仑势Z/r vl ( r )  vloc ( r )   vl ( r ) vˆ PS  vloc (r)  | Ylm  vl (r)Ylm | lm  Kleinman-Bylander可分离格式 KB KB KB | Y   v ( r )  Y |  |   E    lm 1 lm  lm l lm | lm lm 积分数目由约N(N+1)/2个变为N个  ghost state  超软赝势(USPP)  非局域赝势 | i  i | vNL   i | ips  | i   (i  T  vloc ) | ips  O 2p  多能量参考点 Bij  ips |  j  | i    ( B 1 ) ji |  j  j vNL   Bij | i   j | ij 容易验证 (T  vloc  vNL ) | ips   i | ips   模守恒条件成立时，Bij和vNL是厄米的  Qij  i |  j  R  ips |  jps  R  0  除去模守恒条件 PRB 1990, 41, 7892 S  1  Qij | i   j | ij vNL   ( Bij   jQij ) | i   j | ij  | S |  jps  R  i |  j  R ps i  (T  vloc  vNL ) | ips   i S | ips  电荷密度、能量、力需重新定义 投影缀加波方法(PAW)  考虑赝波函数与全电子波函数之间的变换  在原子周围做分波展开  通过求得赝波函数能量,电荷密度等信息 PAW变换  令 ，有 算符变换  对任意局域算符存在对应的赝算符  电荷密度  动能 on site density matrix Hartree能  赝波函数与全电子波函数在球内的模不相等  在处理长程静电相互作用时常引入补偿电荷  Hartree能 单中心赝电荷 单中心补偿电荷 总能  总能可以写成三项之和 赝波函数  赝波函数是PAW方法的变分变量，可通过如下方程自 洽求解  如果分波函数在PAW球内形成完备基组，全电子波函 数与芯态正交，这正是PAW方法被认为是一种全电子 方法的原因 Al 3p 3s 2p 2s 1s effective Al PAW Al 2p 1s 3p 3s nodeless with nodal structure 自洽场迭代  两圈循环优化波 函数和密度  内层循环优化波 函数：CG, DIIS, Davidson算法  外层循环优化电 荷密度：DIIS算 法 多维优化  最速下降法 F ( x ) g |x  xi x xi 1  xi   i1gi  共轭梯度法 f 1 2 f f ( x )  f ( P)   xi   xi x j  ... 2 i , j xi x j i xi 1  c  b  x  x  A x 2 f  Ax  b;  (f )  A x u  [ (f )]  0  u  A  v  0  拟牛顿法 1 f (x  δ)  f (x)  f (x)δ  δAδ 2 f (x  δ)  0  δ   A1f (x) lim Hi  A1 Broyden方法，BFGS 矩阵本征值问题的迭代算法  Krylov子空间方法(Lanczos迭代) n1  cn1 (hˆn  hnnn  hn,n1n1 )  残量 R[ ap ]  (hˆ   ap Iˆ) |  ap  hij  i | hˆ |  j  ap ˆ ap   | h |   ap    ap |  ap   预处理(preconditioning)   (hˆ   ap Iˆ) 1 R[ ap ]  KR[ ap ]  Teter预处理矩阵：对角矩阵，对大的G矢收敛到动能，对 小的G矢趋于一个常量 单带优化  最速下降法  沿残量方向优化本征值等价于一个2x2的本征值问题  每个迭代步引入一个新残量，逐步扩大迭代子空间  共轭梯度法  每步引入的新矢量与之前搜索方向共轭 多带优化  Davidson迭代 同时考虑n个本征矢，引入n个残量  解2n×2n本征值问题，取能量最低的n个波函数   块Davidson方法  在所有带中取一个子集，进行Davidson迭代  所有子集收敛以后，进行一个整体的Raighly Ritz优化 残量优化方法  RMM-DIIS 不优化本征值，而是优化  R[n1 ] | R[n1 ]，  同样对应一个迭代子空间的本征值问题  因为残量的模对每个本征态都是极小值，所以不需要先对 残量做正交归一化  趋向于找到与最初设置相近的极值点，在初值不好的时候 会丢掉一些值  布里渊区积分  电荷密度 occ ＝  n*k  r nk  r  drdk n BZ  不失一般性，积分可以离散成加权求和 1    BZ BZ w k f k  特殊点的选取 f 2 (k )  A0  A1 sin(k )  A2 sin(2k ) I  2 0 f 2 (k )dk  A0  3  f 2 (k  )  f 2 (k  ) 2 2 k Monkhorst-Pack取样法  将三维问题转化为三个一维问题  特殊k点的分数坐标 ur  (2r  nk  1) / 2nk r  1, 2,3,..., nk  对偶数分割不包括原点，以减少计算量  对六角晶格，应该将Gamma点包含在k网格中 金属与绝缘体  占据数零温下为阶跃函数 1   n (k) | Xˆ | n (k) ( nk  )dk n BZ BZ  与绝缘体比较，金属需要更密的k点取样 分数占据(SMEAR)  Fermi-Dirac分布 1 f ( )     e k BT  1  能量对分数占据f不再满足变分原理  新的变分函数是自由能F，包括熵的贡献 F  E  k S ( f nk ) nk Gaussian Smearing  将能级展宽成高斯函数，积分可得 f ( )       1 1  erf    2  k BT    熵不能通过f表达，σ没有物理意义  力通过自由能微分得到，不同于E(0)的力，解决办法： 1 E (T  0)  E0  ( F  E ) 2 Methfessel & Paxton Smearing  将阶梯函数用一组正交函数展开，Gauss是其零阶近似  能量外推  高阶MP方法中F里的熵贡献较小，计算的力更准。  不适用与半导体/绝缘体，占据数可能大于1 Tetrahedron Integration  在由四个k点组成四面体内做线性插 值  对金属误差抵消不完全，调整BZ积分 时k点的权重 (Blöchl Correction)  对分数占据不满足变分原理  best k-point convergence for energy  k meshes for tetrahedra BZ-integration have to include Gamma and the kpoints at the BZ edges Mixing  Linear Mixing niin1   niout  (1   )niin  niin   (niout  niin ) 对束缚很强的刚性系统，可以取较大的 ，而对于金属表面 之类的较软的系统，则比较困难  平面波情形，charge sloshing problem   vout  G    G, G  vin  G   1 G   Kerker Mixing：对介电函数做Thomas-Fermi屏蔽近似 G G G 2 2 G  Gmax 2 min Broyden Mixing  linear mixing可以类比于Quasi-Newton-Raphson迭代 xi+1=xi-J-1Ri  采用J0=αI作为初始的Jacobian，通过quasi-Newton弛豫， 逐步更新J Pulay Mixing  RMM-DIIS like 求几个密度的线性组合系数  Minimize the normal of the residue vector subject to the constraint of conserving the number of electron (DIIS like)  If SCF Convergence Fails 能量泛函  KS能量泛函 EKS  TS [n]  Epot [n] E pot   drvext (r )n(r )  EH [n ]  E XC [n ] TS [n ]  ES   drv (r )n (r ) in out  有效势泛函 EKS [v in ]  ES [v in ]   drv in (r )n out (r )  E pot [n out ]  Harris-Foulkes泛函 EHF [nin ]  ES [vn ]   drvn (r)nin (r)  E pot [nin ] in  in 主要用在非自洽计算时的能量估计，无需计算 nout N ES    i i 几何优化  数学上等价于寻找函数最小值  梯度  牛顿法 最速下降法 the highest frequency mode determines the maximum stable step-width, but the soft modes converge slowest  Number of iterations  DIIS  gradient is linear in it’s arguments for a quadratic function Number of interations  POTIM, NFREE  Trouble Cases for DIIS CG line minisations is done using a variant of Brent algorithm  POTIM  Damped MD  equation of motion  α(POTIM) must be as large as possible, but without leading to divergence  SMASS must be set to (like SD) 几何优化算法选择 上机实践  石墨烯能量收敛测试 K sampling  Ecut  Real space grid   测试mixing参数对石墨烯SCF收敛的影响