密度泛函理论与从头计算分子动力学.pdf
《密度泛函理论与从头计算分子动力学.pdf》由会员分享,可在线阅读,更多相关《密度泛函理论与从头计算分子动力学.pdf(40页珍藏版)》请在三一文库上搜索。
1、密密密度度度泛泛泛函函函理理理论论论与与与从从从头头头计计计算算算分分分子子子动动动力力力学学学 $ 1 引言 自从上世纪60年代以来,密度泛函理论(DFT)建立.并在局域密度近似(LDA)下导出著名 的Kohn-Sham(KS)方程以来,DFT一直是凝聚态物理领域计算电子结构及其特性的最有力 的工具. 近几年来DFT与分子动力学相结合,在材料设计,合成,模拟计算和评价诸多方面有 明显进展,成为计算材料科学的重要基础和核心技术. 特别在量子化学计算领域,1987年以 前主要用Hartree-Fock(HF)方法.但近年来,用DFT的工作以指数增加.以致于HF方法应用已 经相当少.W.Kohn因
2、提出DFT获得1998年诺贝尔化学奖. 已经表明了DFT在计算化学领域 的核心作用与应用的广泛性. 分子动力学计算机模拟是研究复杂的凝聚态系统的有力工具.这一技术既能得到原 子的运动轨迹,还能像做实验一样作各种观察.对于平衡系统,可以在一个分子动力学观 察(observation time)内作时间平均来计算一个物理量的统计平均值.对一个非平衡系统过 程,只要发生在一个分子动力学观察时间内(1-10ps)的物理现象也可以用分子动力学计算进 行直接模拟.可见数值实验对理论与实验的有力补充,特别是许多与原子有关的微观细节,在 实验中无法获得而在计算机模拟中可以方便的得到. DFT与分子动力学(MD
3、)相结合可以有大量不同类型的应用. 如晶格生长,外延生长,离 子移植,缺陷运动,无序结构,表面与界面的重构,电离势的计算,振动谱的研究,化学反应的 问题,生物分子的结构,催化活性位置的特性以及材料电子结构和几何结构.固,液体的相变 等.现在这些方法已发展成为成熟的计算方法.DFT 的另一个特点是,它提供了第一性原 理(fi rst-principal)或称为从头计算的理论框架.在这个框架下可以发展各种各样的能带计 算方法.虽然在DFT的所有这些实际应用中,几乎都采用局域密度近似(LDA),这是一种不能 控制精度的近似,因而DFT方法的有效性在很大程度上要看结果与实验的一致性.人们没 有任何直接
4、的方法可以改善LDA的精度. 然而DFT允许发展别的方法加以补充.例如广义 梯度近似(GGA)等方法,把密度分布n( r )的空间变化包括在方法之中,实现了比较大幅度减 少LDA误差的目的. 相比较传统的量子化学方法,如组态相互作用(CI)方法,DFT+MD方法显然可以用于数 百个原子的大分子体系.但对于具有强关联相互作用的体系,似乎化学家更愿意应用CI方 法.而对于凝聚态物理领域,DFT+MD方法可以相当精确的计算材料的电子结构和相应的 许多物性. 在DFT获得巨大成功的同时,也有些不可忽视的弱点与困难.针对这些问题已经发展了许 2 多不同的方法,这些方法可以用Kohn-Sham方程的有效H
5、amiltonian的各部分和波函数构造 上的考虑进行归类,如图一所示 传统的DFTLDA是预言多电子体系基态性质的理论。对于激发态性质的描述总是与 实验不符合,这固然因为(1)激发态本身存在着复杂性质。(2)DFTLDA理论存在 着对激发态描述困难。此外,对于具有强关联相互作用体系LDA理论描述不够好。因而 发展出LDA+U和LDA+等方法。此外关于大原子数的复杂体系,近几年来发展了各种 线性标度的方法,也称为O(N)算法。为研究复杂体系提供了有力的工具。 $ 2 基本理论 首先,对包括诸多电子,离子,原子的凝聚态体系来说,面临的首要问题是,由于原 子核较重而且运动中无法跟随电子的运动。因而
6、可以把电子与原子分开考虑,原子实作为 3 带有正电荷的外场存在。这就是所谓的绝热近似或称BO近似(BornOppenheimer).因而 我们可以把电子的哈密顿量写成 H = T + V + Vext(2.1) $ 2.1 Hohenberg-Kohn定理 1964年,HK提出HK两个定理: 1) 多电子系统在外场势Vext作用下,其基态的电子密度( r), 与基态下任意力学量O的 可观测量是一一对应的,且是严格的基态电子密度( r)的泛函. =O(2.2) 2) 若O为哈密顿量H,则基态的总能量函数H=EV ext具有如下形式 EV ext=+(2.3) 而对任意多电子体系的普适量FHK为
7、FHK= EV ext=FHK+R( r)Vext( r)d r(2.4) 这里,我们不准备证明上述两个定理,只说明如下3点 1)任何与力学量或势之间存在一一对应性(one to one correspondence) 2)普适函数FHK可以很容易用密度算符表示。 3)HK第二定理可以使我们利用变分定理,通过求能量极小值来获得基态的电子密度。 $ 2.2 Kohn-Sham方程 KS方程诞生于1965年,KS方程提出后,才真正使DFT成为实际计算可能。设总能 量函数Ee (包括电子间相关能)及在Hartree-Fock近似下的总能量函数为EHF(不包 括电子间相关能) 则Ee=T+V(2.5)
8、 EHF= T0+ (VH+ VX |z V ) (2.6) T,V为严格的动能及电子间的势能函数。T0为无相互作用的电子气动能,VH为电子间 的库仑能,VX为电子间交换能。因而,不难得到上述两种体系下能量差别仅仅在于电子 相关能VC. 则 VC=Ee EHF=T-T0(2.7) 则HK的普适函数FHK可以写为: FHK= T+V +T0T0= T0+V +(T T0 |z VC ) = T0+V +VC+VHVH= T0+VH+VC+(V VH |z VX ) 4 =T0+VH+(VX+ VC |z VXC )(2.8) 这里VXC被称为交换关联势函数,这样总能量泛函可以写为 EVext =
9、 T00 + (1/2) Z d rdr0 o( r)o(r0) | r r0| + Z Vextd r + EXCo (2.9) 这里交换关联势可由下式给出 VXC= EXC 对应的系统哈密顿量(通常称为KS哈密顿量)写为 Hksi= 1 2 2 i + Veffi= Eii(2.10) 这里单粒子波函数满足 ( r) = N X i=1 i( r)i( r) (2.11) Veff= Vext( r) + Z d r (r0) | r r0| + EXC ( r) (2.12) (2.9)(2.12)式就构成了KS方法的基本方程组。 遗憾的是,它还是形式上的东西,因为包含电子多体效应交换关
10、联效应的EXC的具体 表达式还不清楚,需作某种近似.现在广泛应用的是局域密度近似(LDA),即把非均匀 电子系统分割成一些小块。在这些小块中认为电子分布是均匀的。这样 r处的子块中的交 换关联能密度xc( r)只取决于该点的( r) ,整个系统的交换关联能为 EXC= Z d r( r)xc( r)(2.13) 至 于 交 换 关 联 势EXC的 具 体 表 达 式 不 断 有 著 者 给 出 。 目 前 用 的 最 广 泛 的 是Cepeley和Alder在八十年代用Monte Carlo方法导出的解析表示(Phys.Rev.lett(566- 568)Vol.45.No.7)。从方法推导近
11、似来看,LDA只局限于那些电子密度缓变的系统。因 而,在1996年,Perdew,Burke,Ernzerhof把它推广到包括电子梯度密度在内的函数,并给 出了解析表示式。称为GGA方法。(Generalized Gradient Approximation.Phys.Rev.lett Vol.77.No.8(3895-3897),如包含电子自旋极化(LSD)情况在内,则(2.13)可写为 ELSD XC n ,n = Z d3r n unif xc n ,n (2.14) 而在GGA下,则为 ELSD XC n ,n = Z d3r f(n ,n ,n ,n )(2.15) 5 现在GGA方
12、法已经得到广泛的应用。总之,局域密度泛函理论,自上世纪60年代提出以 来,经过许多人的发展与完善,已广泛的应用于各领域多原子体系的电子结构计算。由于 在此理论近似下,单电子运动方程中的交换关联势形式简单,使计算工作量大大减少。因 此,现在大多数能带或者原子集团模型中的计算方法,都是在局域密度泛函理论近似的基 础上建立起来的。 $ 2.3 Kohn-Sham方程的求解 在传统的固体,分子的局域密度泛函电子结构计算中,KS自洽方程组一般可以这样求 解:先给定一个电荷o( r),比如说可以选用各原子电荷密度的迭加,然后代入(2.12)计算有 效势Veff,再通过对角化来求解(2.10) ,得到HKS
13、的本征矢。即i( r)。然后可以求得新的分 布电荷密度o( r),继续进行迭代直致自洽。由于对角化计算工作量以M3增加(M为KS轨 道展开所需基函数的个数),因而对有数百个原子的系统工作量非常庞大。 KS方程组也可以通过其它方法解。为了求得能量泛函对电子密度变分的最小值,我们 可以在电子自由度i的空间中引入一个虚拟的动力学过程,来使泛函达到极小值。最简单 动力学考虑是速降法。 i( r,t) = E i( r,t) (2.16) 上式中的点 代表对时间的微分。变分 E i = HKSi,依赖于i,求解上式时,为保证波 函数的正交性,还得加上一个约束条件,虚拟动力学方程是: i( r,t) =
14、E i( r,t) + X j ijj(2.17) 实际上,先用无约束的动力学方程(2.16)计算i,然后通过适当的正交化过程对i进行 修正,并得到系数ij。当系统达到极小值时,i= 0 ,即HKSi= P jijj,把相应的约束 矩ij对角化后,既可以得到KS方程的本征矢与本征值。经验表明,对于一定的体系坐 标E只有一个极小值。这样速降法也能得到能量极小值。这种方法与传统的方法相比省去 了大量的对角化计算。 速降法的有效性取决于波函数达到收敛的有效步数,这个步数可能很多,特别是在低 对称下,人们也做了许多改善速降法的尝试,如有更复杂的积分方程,以及改善的速降 法,及共扼梯度法等等。 象速降法
15、及共扼梯度法这样的算法在哈密顿量HKS作用到i的每步计算中都需要适当 的正交化操作。如果i是用平面波展开的(设N个单电子轨道均用M个平面波展开), 则解HKSi就需要NM2次浮点计算。加上正交化需要N2M次操作,计算量很大。为了减少 计算量,可以把算符HKS分成动能项与势能项,动能项在倒格矢空间是对角的,作用到 6 波函数i上只需O(NM)次操作。势能Veff是局域的,则作用到波函数i上也可方便的用 快速付立叶变换技术(FFT)在实空间计算, 也只要O(NMlnM)次操作。因而工作量以线性 或O(NMlnM)次增长,这与传统的标准对角化技术相比,工作量有显著减小. $ 2.4 局域密度近似下的
16、离子间相互作用势的计算 传统的,如果给定原子坐标 RI,由从头计算方法推导离子间相互作用势V RI,然后 进行第一性原理的分子动力学模拟,一般有三个基本步骤,即在每个MD 过程中:(1)对于给定 的 RI求解自洽的KS方程.(2)根据HellmanFeymam的力定理计算每个离子受到的力.(3) 求解牛顿运动方程 MI RI= RIV 但需要同时对电子和离子的自由度进行优化(Bendt,Zunger).利用速降法,通过下面的方 程组就可以简单的实现这样的想法. ii= E i + X j ijj(2.18a) MI RI= E RI (2.18b) 其中约化“质量”i和Mi的引入可以调节与两套
17、参数i和 RI相关的时间标度. 方 程组(2.18)就可以确定从能量面上初始点对附近的极小值之间的一条轨迹. 在 RI不变 时,通过变化i,就可以找到能量面上E的唯一的极小值. 相应的在BO势能面上也有与之 对应的极小值;V RI=min |z i Ei, RI , 这样能量面E上的局域极小值也就成了势能 面V上的局域极小值.通过同时对i, RI进行变分,方程组(2.18)就可以用于复杂的原子 结构的计算.与传统的第一性原理方法直接应用于分子动力学模拟相比,这一方法在计算量 上大为减少. 不过,方程组(2.18)决定的动力学过程还有很大的局限性.它所得到的结果只是局域优化 的结果.实际上, 势
18、能面V一般有许多个极小值,为了解决这个问题,1985年,Car,Parrinello提 出一个比较满意的全局优化算法. $ 3 从头计算分子动力学(CP方法) $ 3.1 CP方法原理 Car,Parrinello提出的从头计算分子动力学方法,最重要的一点是在真实的物理系统中 引入一个虚拟的电子动力学系统,这样组成的新系统的势能面E是离子与电子的一个总泛 函.使得这个虚拟系统产生的轨迹与具有同样势能面V的实际物理系统的轨迹一致.并能同 时处理电子离子系统. 7 实际物理系统的经典拉氏量为: LCL= 1 2 P IMI R2 I V RI(对离子晶格而言)(3.1) 由拉氏量L=T-V,则可以
19、把虚拟系统的广义经典拉氏量写为: L = occ X i Z d ri| i|2+ 1 2 X I MI RI 2 ERI,i + X i,j ij Z d r( ij ij) (3.2) 这里,L为两套自由度i, RI的泛函,它本身虽不显含时t,但i,RI是均与时间有关.i为可 调参数.(单位:质量 长度长度),相当于电子的“质量”实际上起着调节电子时间运动标 度的作用. 可令i= 与波函数无关.其中Ei, RI应为电子和离子耦合的能量函数,应 包括电子动能,电子相互作用能(Hartree能),电子的交换关联能. 以及离子实对电子的作用 能,以及离子间的相互作用能.即等价于在计算中的所需的“
20、总能量”(即除了离子动能 项之外)。拉格朗日乘子ij是为了保证i的正交性而引入的。在经典力学中,就是一个 完整约束。根据拉格朗日-欧拉方程: L qi d dt L qi = 0i = 1.2.3n. 从(3.2)式不难得到: i= E i + X j ijj(3.3a) MI RI= E RI (3.3b) 我们可以通过改变离子“速度” RI和电子“速度” i,可以控制这一虚拟系统的温度。 从而可以对它们进行各种热力学准静态过程处理,如退热,冷却等。在这一热力学过程 中,可以同时处理离子和电子的自由度。特别可以设计一个缓慢“降温”的过程,使系 统T 0K时,达到平衡状态,这一方法被称为动力学
21、退热模拟。 用动力学退热模拟来做全局优化计算的优点十分明显.在有限温度下,各自由度的速度是 不为零的, 而且由于系统的复杂性,会产生随机的热力学运动.由各态经历可知,只要初始温 度足够高,体系就能达到势能面上的各个点. 通过降温,体系在势能面上随机运动慢慢的局 限于几个极小值之间的跳跃.而且跳跃的频率逐渐减小.只要温度冷却的足够慢, 系统最终 一定会平衡在势能面的最低点.但有时这种全局优化方案在系统即将达到能量最低点时, 其 冷却速度会慢的无法忍受.以至于在实际操作中,计算机无法收敛,此时,就需要其他因素来 进行调节. 上述我们提到的温度只是离子的“温度”,并没有涉及与电子相关的温度.在整个动
22、力学 模拟过程中,电子自由度与离子自由度之间并不是始终处于热力学平衡状态,还可能会使按 8 这种方法模拟的离子运动轨迹与实际物理系统中离子的运动轨迹不一致. 离子的实际运动 方程为: MI RI= V RI RI (3.4) 由于V RI=min |z i Ei, RI 相关,所以要使(3.3b) 与(3.4) 产生离子轨道一致,要 求在模拟过程中,泛函ERI,i 始终处于对i的极小值.我们可以通过参数 和初始条 件i ,i的选取,使电子运动时间标度远小于离子运动时间标度(BO条件) ,以使电子在 离子坐标每次变动之前尽量趋于基态,结果电子自由度在能量面E的极小值附近做简谐运 动,从而使离子的
23、轨迹在高温和低温下基本上都保持在BO面上.即若和i ,i的选取使 离子和电子的自由耦合度很弱, 它们之间能量转移是较小,而使电子的运动绝热于离子的运 动,电子可以调整自己的运动状态来“追随”离子的运动. 这样除了对初始的原子构型外,不 必每一步自洽求解KS方程组,从而大大减少了工作量.真正使第一性原理用于分子动力学模 拟. 现在,这一方法已有效的运用到如Cluster,液态金属, 无序系统等许多方面. $ 3.2 基矢与原胞 最初CP提出的从头计算分子动力学方法是在赝势和平面波的基础上具体实现的.(当然 现在也有用缀加平面波(augmented plane wave)的算法.如软件包WIN2K
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 密度 理论 从头 计算 分子 动力学
链接地址:https://www.31doc.com/p-3705673.html