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