密度泛函理论与从头计算分子动力学.pdf

密度泛函理论与从头计算分子动力学$1引言自从上世纪60年代以来,密度泛函理论(DFT)建立.并在局域密度近似(LDA)下导出著名 的Ko hn-Sham(KS)方程以来,DFT一直是凝聚态物理领域计算电子结构及其特性的最有力 的工具.近几年来DFT与分子动力学相结合，在材料设计，合成，模拟计算和评价诸多方面有 明显进展,成为计算材料科学的重要基础和核心技术.特别在量子化学计算领域1987年以 前主要用Har t r ee-Fo c k(HF)方法.但近年来用DFT的工作以指数增加.以致于HF方法应用已 经相当少.W.Ko hn因提出DFT获得1998年诺贝尔化学奖.已经表明了DFT在计算化学领域 的核心作用与应用的广泛性.分子动力学计算机模拟是研究复杂的凝聚态系统的有力工具.这一技术既能得到原 子的运动轨迹，还能像做实验一样作各种观察.对于平衡系统何以在一个分子动力学观 察(o bs er v at io n t ime)内作时间平均来计算一个物理量的统计平均值.对一个非平衡系统过 程，只要发生在一个分子动力学观察时间内(Ll Ops)的物理现象也可以用分子动力学计算进 行直接模拟.可见数值实验对理论与实验的有力补充，特别是许多与原子有关的微观细节，在 实验中无法获得而在计算机模拟中可以方便的得到.DFT与分子动力学(MD)相结合可以有大量不同类型的应用.如品格生长孙延生长，离 子移植，缺陷运动，无序结构，表而与界面的重构，电离势的计算，振动谱的研究化学反应的 问题，生物分子的结构，催化活性位置的特性以及材料电子结构和几何结构.固，液体的相变 等.现在这些方法已发展成为成熟的计算方法.DFT的另一个特点是，它提供了第一性原 理(fir s t-pr inc ipal)或称为从头计算的理论框架在这个框架下可以发展各种各样的能带计 算方法.虽然在DFT的所有这些实际应用中，几乎都采用局域密度近似(LDA),这是一种不能 控制精度的近似，因而DFT方法的有效性在很大程度上要看结果与实验的一致性.人们没 有任何直接的方法可以改善LDA的精度.然而DFT允许发展别的方法加以补充.例如广义 梯度近似(GGA)等方法，把密度分布n)的空间变化包括在方法之中，实现了比较大幅度减 少LDA误差的目的.相比较传统的量子化学方法，如组态相互作用(CI)方法QFT+MD方法显然可以用于数 百个原子的大分子体系但对于具有强关联相互作用的体系，似乎化学家更愿意应用CI方 法.而对于凝聚态物理领域QFT+MD方法可以相当精确的计算材料的电子结构和相应的 许多物性.在DFT获得巨大成功的同时，也有些不可忽视的弱点与困难.针对这些问题已经发展了许2Self-consistent All-electron full potential Harris functional-Ail-electi-cn mufnn-tin-Pseudopotential-JelliumGW.-Time-dependent DFT Beyond LDA-Optimized effective potentials(OEP)LDA+U,LDA+一 Self-interacticn correcions(SIC)Fully-relativistic-Semi-relativistic-Non-relativistic-Local density approximation(LDA)4v2+K(r)+(r)iirf 二E：必Non-periodic-Periodic-Symmetry-Locality,O(N)一 Non-spin-polarized Spin-polarized-Atomic orbitals Slater type(STO)I Gaussians(GTO)*-Numerical(Dmol,TB)-1-Plane waves(FLAPW)Spherical waves(LMTO,ASW)-Fully numerical-O(N)一吨Plane waves-Augmentation图1密度泛函理论的发展多不同的方法，这些方法可以用Ko hn-Sham方程的有效Hamil t o nian的各部分和波函数构造 上的考虑进行归类,如图一所示传统的DFT LDA是预言多电子体系基态性质的理论。对于激发态性质的描述总是与 实验不符合，这固然因为（1）激发态本身存在着复杂性质。（2）DFT LDA理论存在 着对激发态描述困难。此外，对于具有强关联相互作用体系LDA理论描述不够好。因而 发展出LDA+U和LDA+等方法。此外关于大原子数的复杂体系，近几年来发展了各种 线性标度的方法，也称为0（N）算法。为研究复杂体系提供了有力的工具。$2基本理论首先，对包括诸多电子，离子，原子的凝聚态体系来说，面临的首要问题是，由于原 子核较重而且运动中无法跟随电子的运动。因而可以把电子与原子分开考虑，原子实作为 3带有正电荷的外场存在。这就是所谓的绝热近似或称B0近似(Bo r n-Opp6nheimer).因而 我们可以把电子的哈密顿量写成H=T+V+Vext(2.1)$2.1 Hohenberg-Kohn定理1964年，HK提出HK两个定理：1)多电子系统在外场势匕J乍用下，其基态的电子密度p(f%与基态下任意力学量。的 可观测量是一一对应的，且是严格的基态电子密度p(竹的泛函.=0p(2.2)2)若。为哈密顿量方，则基态的总能量函数Hp=Ev p具有如下形式EveXtp=+(2.3)而对任意多电子体系的普适量尸HK，为FHKp=EVext p=Fhk p+fp vext(r)dr(2.4)这里，我们不准备证明上述两个定理，只说明如下3点1)任何p与力学量或势之间存在一一对应性(o ne t o o ne c o r r es po ndenc e)2)普适函数尸可以很容易用密度算符表示。3)HK第二定理可以使我们利用变分定理，通过求能量极小值来获得基态的电子密度。$2.2 Kohn-Sham方程K S方程诞生于1965年,K S方程提出后，才真正使DFT成为实际计算可能。设总能 量函数乙同(包括电子间相关能)及在Har t r ee-Fo c k近似下的总能量函数为石印加(不包 括电子间相关能)贝(J Ee=T+V(2.5)Ehf=Tq+(VH+Vx)(2.6)VT,V为严格的动能及电子间的势能函数。毁为无相互作用的电子气动能，质为电子间 的库仑能，口为电子间交换能。因而，不难得到上述两种体系下能量差别仅仅在于电子 相关能.则Vc=Ee-Ehf=T-Tg(2.7)则HK的普适函数外长可以写为：Fhk T+V+TqTq Tq+V+(T 7q)Tq+V-Vc-Vh Vh Tq+Vh+Vc+(V Vn)Vc Vx4=n+0+(“+%)(2.8)-一 Vxc 这里“被称为交换-关联势函数，这样总能量泛函可以写为EVextp=T0p0+(1/2)/而加。言索+/Vextpdr+Exc p。(2这里交换关联势可由下式给出SEXcp yxc=7-而对应的系统哈密顿量(通常称为KS哈密顿量)写为ks6i=-|v +=E(2.10)这里单粒子波函数满足N，(了)=。:&(2.11)2=1Veff=Vextf+/*产2T+笔F(2.12)J f-rf M(/)(2.9)-(2.12)式就构成了KS方法的基本方程组。遗憾的是，它还是形式上的东西，因为包含电子多体效应交换关联效应的Exc的具体 表达式还不清楚，需作某种近似.现在广泛应用的是局域密度近似(LDA),即把非均匀 电子系统分割成一些小块。在这些小块中认为电子分布是均匀的。这样F处的子块中的交 换关联能密度”只取决于该点的p/)，整个系统的交换关联能为Exc=/历)5p(2.13)至于交换关联势Exc的具体表达式不断有著者给出。目前用的最广泛的 是Cepel ey和Al der在八十年代用Mo nt e Car l o方法导出的解析表(P hys.R ev.l et t(566-568)Vo l.45.No.7)o从方法推导近似来看，LDA只局限于那些电子密度缓变的系统。因 而，在1996年，P er dew,Bur kgEr nzer ho f把它推广到包括电子梯度密度在内的函数，并给 出了解析表示式。称为GGA方法。(Gener al ized Gr adient Appr o ximat io n.P hys.R ev.l et t Vo l.77.No.8(3895-3897),如包含电子自旋极化(LSD)情况在内，贝(J(2.13)可写为E射n。八=/d3r n嘤订n=八(2.14)而在GGA下，则为E肥 n?,ni=I d3r f(n n Vn Vn(2.15)5现在GGA方法已经得到广泛的应用。总之，局域密度泛函理论，自上世纪60年代提出以 来，经过许多人的发展与完善，已广泛的应用于各领域多原子体系的电子结构计算。由于 在此理论近似下，单电子运动方程中的交换关联势形式简单，使计算工作量大大减少。因 此，现在大多数能带或者原子集团模型中的计算方法，都是在局域密度泛函理论近似的基 础上建立起来的。$2.3 Kohn-Sham方程的求解在传统的固体，分子的局域密度泛函电子结构计算中，KS自洽方程组一般可以这样求 解：先给定一个电荷p/竹,比如说可以选用各原子电荷密度的迭加，然后代入(2.12)计算有 效势匕万,再通过对角化来求解(2.10),得到方ks的本征矢。即&。然后可以求得新的分 布电荷密度外(可，继续进行迭代直致自洽。由于对角化计算工作量以河3增加(M为KS轨 道展开所需基函数的个数)，因而对有数百个原子的系统工作量非常庞大。KS方程组也可以通过其它方法解。为了求得能量泛函对电子密度变分的最小值，我们可以在电子自由度出的空间中引入一个虚拟的动力学过程，来使泛函达到极小值。最简单 动力学考虑是速降法。(2.16)上式中的点代表对时间的微分。变分患=HksM依赖于仇，求解上式时，为保证波 函数的正交性，还得加上一个约束条件，虚拟动力学方程是：6”二+工 血 17)实际上，先用无约束的动力学方程(2.16)计算所,然后通过适当的正交化过程对仇进行 修正，并得到系数入户当系统达到极小值时，加=0,即为5弧=2他与把相应的约束 矩他对角化后，既可以得到KS方程的本征矢与本征值。经验表明，对于一定的体系坐 标E只有一个极小值。这样速降法也能得到能量极小值。这种方法与传统的方法相比省去 了大量的对角化计算。速降法的有效性取决于波函数达到收敛的有效步数，这个步数可能很多，特别是在低 对称下，人们也做了许多改善速降法的尝试，如有更复杂的积分方程，以及改善的速降 法，及共扼梯度法等等。象速降法及共扼梯度法这样的算法在哈密顿量方ks作用到&的每步计算中都需要适当 的正交化操作。如果仇是用平而波展开的(设N个单电子轨道均用M个平面波展开)，则解就需要N2次浮点计算。加上正交化需要N2M次操作，计算量很大。为了减少 计算量，可以把算符Hks分成动能项与势能项，动能项在倒格矢空间是对角的，作用到 6波函数而上只需O(NM)次操作。势能匕万是局域的，则作用到波函数仇上也可方便的用 快速付立叶变换技术(FFT)在实空间计算，也只要O(NMl nM)次操作。因而工作量以线性 或O(NMl n 次增长，这与传统的标准对角化技术相比，工作量有显著减小.$2.4局域密度近似下的离子间相互作用势的计算传统的，如果给定原子坐标%，由从头计算方法推导离子间相互作用势V必,然后 进行第一性原理的分子动力学模拟，一般有三个基本步骤，即在每个MD过程中:对于给定 的%求解自洽的KS方程.根据Hel l man-Feymam的力定理计算每个离子受到的力.求解牛顿运动方程=-7RlV但需要同时对电子和离子的自由度进行优化(Bendt,Zunger).利用速降法,通过下面的方 程组就可以简单的实现这样的想法.口加=一菽+E 血(2.18q).3 E=-(2.186)oRj其中约化“质量”小和方z的引入可以调节与两套参数6和 R 相关的时间标度.方 程组(2.18)就可以确定从能量而上初始点对附近的极小值之间的一条轨迹.在 R/不变 时，通过变化&,就可以找到能量而上E的唯一的极小值.相应的在BO势能面上也有与之 对应的极小值;V Ez=mm Ea,R,这样能量而E上的局域极小值也就成了势能 皿面V上的局域极小值.通过同时对&,R 进行变分,方程组(2.18)就可以用于复杂的原子 结构的计算.与传统的第一性原理方法直接应用于分子动力学模拟相比，这一方法在计算量 上大为减少.不过,方程组(2.18)决定的动力学过程还有很大的局限性.它所得到的结果只是局域优化 的结果.实际上，势能而V一般有许多个极小值,为了解决这个问题1985年。ar Far r inel l o提 出一个比较满意的全局优化算法.$3从头计算分子动力学(CP方法)$3.1CP方法原理Car far r inel l o提出的从头计算分子动力学方法，最重要的一点是在真实的物理系统中 引入一个虚拟的电子动力学系统,这样组成的新系统的势能面E是离子与电子的一个总泛 函.使得这个虚拟系统产生的轨迹与具有同样势能面V的实际物理系统的轨迹一致.并能同 时处理电子-离子系统.7实际物理系统的经典拉氏量为：LCL=|/Mi瑶-（对离子品格而言）（3.1）由拉氏量1=则可以把虚拟系统的广义经典拉氏量写为：(3.2)这里,L为两套自由度&,必的泛函，它本身虽不显含时t,但&园是均与时间有关心为可 调参数.（单位:质量x长度x长度），相当于电子的“质量”实际上起着调节电子时间运动标 度的作用.可令的=4与波函数无关.其中E 6,%应为电子和离子耦合的能量函数，应 包括电子动能，电子相互作用能（Har t r ee能）,电子的交换关联能.以及离子实对电子的作用 能，以及离子间的相互作用能.即等价于在计算中的所需的“总能量”（即除了离子动能 项之外）。拉格朗日乘子八刀是为了保证仇的正交性而引入的。在经典力学中，就是一个 完整约束。根据拉格朗日-欧拉方程：i=1.2.3-n.从（3.2）式不难得到：b i=菽+52 人力功（3.3）.dEMiRi=一万元-（3.36）Orii我们可以通过改变离子“速度”A/和电子“速度”而，可以控制这一虚拟系统的温度。从而可以对它们进行各种热力学准静态过程处理，如退热，冷却等。在这一热力学过程 中，可以同时处理离子和电子的自由度。特别可以设计一个缓慢“降温”的过程，使系 统T-OK时，达到平衡状态，这一方法被称为动力学退热模拟。用动力学退热模拟来做全局优化计算的优点十分明显.在有限温度下，各自由度的速度是 不为零的，而且由于系统的复杂性，会产生随机的热力学运动.由各态经历可知，只要初始温 度足够高,体系就能达到势能面上的各个点.通过降温,体系在势能面上随机运动慢慢的局 限于几个极小值之间的跳跃.而且跳跃的频率逐渐减小.只要温度冷却的足够慢，系统最终 一定会平衡在势能而的最低点.但有时这种全局优化方案在系统即将达到能量最低点时，其 冷却速度会慢的无法忍受.以至于在实际操作中，计算机无法收敛，此町就需要其他因素来 进行调节.上述我们提到的温度只是离子的“温度”，并没有涉及与电子相关的温度.在整个动力学 模拟过程中，电子自由度与离子自由度之间并不是始终处于热力学平衡状态，还可能会使按8这种方法模拟的离子运动轨迹与实际物理系统中离子的运动轨迹不一致.离子的实际运动 方程为：访的=一驾a(3.4)0R1由于V Ri=mm E6,R,相关，所以要使(3.3b)与(3.4)产生离子轨道一致，要 依求在模拟过程中，泛函始终处于对济的极小值.我们可以通过参数4和初始条 件6，戒的选取，使电子运动时间标度远小于离子运动时间标度(B。条件)，以使电子在 离子坐标每次变动之前尽量趋于基态，结果电子自由度在能量而E的极小值附近做简谐运 动，从而使离子的轨迹在高温和低温下基本上都保持在B0而上.即若从和6，4的选取使 离子和电子的自由耦合度很弱，它们之间能量转移是较小，而使电子的运动绝热于离子的运 动，电子可以调整自己的运动状态来“追随”离子的运动.这样除了对初始的原子构型外不 必每一步自洽求解KS方程组，从而大大减少了工作量.真正使第一性原理用于分子动力学模 拟.现在,这一方法已有效的运用到如Cl us t er,液态金属，无序系统等许多方而.$3.2基矢与原胞最初CP提出的从头计算分子动力学方法是在震势和平面波的基础上具体实现的.(当然 现在也有用缀加平面波(augment ed pl ane w av e)的算法.如软件包WIN2K.但主流是基于度 势和平面波的，主要是因为平面波基能方便的采用FFT技术,使能量，力等的计算在实空间和 倒空间快速变换，这样计算尽可能在方便的空间中进行.如前面讲的哈密顿量中动能项的矩 阵元，在倒空间中只有对角元非零.就比实空间减少了计算量.此外，平而波基函数具体形式 并不依赖于离子坐标，这样，一方而，价电子对离子作用力可以直接应用Hel l man-F6ymann定 理得到解析表示，另一方而也使总能量的计算在不同原子位型下有基本相同的精度.有利 于分子动力学模拟.此外，平而波的计算的收敛性和精确性比较容易控制，这可以通过截断 能(Ecut)来控制平而波基的数量.当然平面波基也有缺点.一般电子轨道都有一定局域性，而 平面波是空间均匀的.因此，电子轨道展开时,与其它原子轨道波函数相比，平而波基的数量 要多得多,这样为了尽量减少平而波基的数量，需要引入鹰势来描写离子实与电子间相互作 用.使电子轨道波函数在离子实内部变化尽量平缓.分子动力学模拟中需要一个满足周期性边界条件的原胞(MDBo x)对于理想晶体的 计算，还是很自然的，因为哈密顿量具有平移对称性，只取其一个单胞即可.但对于无序系 统(无定型结构的固体与液体，表而，界而缺陷及杂质材料等问题，就必须把原胞取的足够 大,Mako v与P ayne已经证明:对于非周期性结构的体系，若用周期性边界条件则能量收敛 9为-5,L为计算中选取的超胞的尺度.因此，这种基于局域密度近似并利用平移对称性来计 算电子结构的方法,对有序及无序都适用，这也是它的优点所在.采用周期性边界条件后，单粒子轨道波函数满足Bl o c h定理，可用平而波展开为.(A)=exp(ifc r)G)exp(zGr)(3.4)G为原胞倒格矢了是简约布里渊内的波矢。(冗份为单粒子轨道波函数的展开系数.(3.4)中 对力的求和可以截断成有限的,给定一个能量生对守的求和可以限制在(l/2)(fe+G)2 2)都很小，可以认为，2肛与/无关,一般就 取相同的形式.NCP P具有如下的性质：1)能量本征值与实际本征值一致.2)当鹿(芯半径)时腐波函数嬴=%(真实单电子波函数).要求度波函数及其一，二级微商在r=晨时连续.3)当底时，震波函数无节点，而总电荷密度与全电子波函数在离子实内分布相同.即 6n6nc=6n6n小4)要求满足 医 一(1/2)(丁。九)7 nr=rc=r=，cCtrjCLT上述第1.2两条是为了满足计算的要求，而3.4两条则为了使NCP P对周围不同化学环境有很 强的普适性.10HSC把离子实的震势丹加分解成两部分，与1无关的匕”e及与1有关的吟6m=Ko r e+/加(3.6)7 2匕。e=-c re erfc ay冏 r i=i(3.7)这里er fc(x)为余误差函数,er f(x)为误差函数,z/5 2erf cx)=1 erf(x)=1-J ez dz(3.8)3Kon=(4+/24+3)6-(3.9)i=l这里Of*芯。？4)4十3回均表示不同原子序数,在NCP P相应的数据库中给出.采 用NCP P后计算所得到的原子震波函数与实际原子波函数相比，的确平坦了许多.如 图2所示.这样,有利用较少数量平面波基的展年便于在实空间和倒空间作快速FFT变换.NCP P自 从提出以来得到广泛应用，但在具有强局域性价轨道的元素的应用中遇到困难.因而才会 有Vander bil t提出的超软度势的诞生.$3.4 超软势(Ul t r as o ft P s eudo po t ent ial USP P)HSC等发展起来的第一性原理范数不变度势,在局域密度近似的基础上用平而波基做电 子结构计算取得很好效果.然而在一些具有强局域性的轨道(主要是元素周期表中第一行元 素和过渡金属元素)体系的应用中遇到困难.例如对。2P轨道，如图3所示，在芯域内部，全电子波函数本身无节点，所以在范数不变(即芯域内的总荷量守恒)条件 下JS波函数不可能比全电子波函数平滑许多.90年代,D.Vander bil t提出的超软鹰势是一种完全非局域的鹰势，也是第一性原理的潢势.本 身具有分离的特点，方便有关能量的计算.超软震势的最大优点在于去掉了原先的范数不变 的约束条件，能使鹰波函数在芯内尽可能平缓，而不致于影响其对各种化学环境的普适性.$341USPP的构造D.Vander bil t第一步先证明完全非局域的鹰势(KB型)，可以通过波函数直接构造.对一个 自由原子用第一性原理全电子计算解电子轨道的波函数和本征值，求解Sc hr o dinger方程小+心石，=&口(3.10)11lines)and a b initio full-core atomic valence wave functions(broken lines)for Mo.The lower panel shows the corresponding bare-ion pseudopotentials.i表示全量子数，而实际上才目应原子数据表上价和.都很容易得到，从而可以 局域密度近似得到全电子的有效势吗石.然后来构造鹰原子波函数，设“为原子实 芯域上的截断半径，则要求我与口在r=r“处光滑连接同时要求满足范数不变条件.R=(2 十4o c)|电 (3.11)12r(a.u.)FIG.l.Oxygen 2p radial wave 吗；嘿,漆黑翼 CSponding pseudo-wave-functions generated using HSC(dot ine)and current(dashed line)methods.同时满足(，R)时此=0,(r 2)时，选定的有效局域势Vloc=心处并且 在口底时，与心石光滑连接.具体选法请参阅P hys,R ev.Vo l.41.No.n(1892)然后,据此可定 义一个非局域的鹰势算符心力Vnl=北%I(3.12)由矩阵元仇儿)可以定义一套局域基函数区=(&%(3.13)3这样心力也可以由股 定义为五乙=&，股=工=工场J XmiXm 宙 ij ij a n13=（召 1）mJ|Xm 1）A坛=2 mjXm ji-Xi ij m n ij从而（T+VLOc+VnlMG=第二步D.Vand6r bi止为了使震波函数具有更好的散射性质（s c at t er ing pr o per t ies）,即使度 波函数对能量的微分在更多的能量本征值3处与全电子波函数一致，加强了范数不变的条 件。在第一步选取鹰波函数的基础上，还要求&满足广义的范数不变条件Q,三0其中Q,=R R（3.15）可以证明，当，j=0时，矩阵场，是厄密的。因而算符VL也是厄密的。从而在构造震势 时，通过增加能量参数本征值.，就可以使鹰势在一定的能量范围内与全电子度势拟合到所 需的精度。设&（）的径向部分为尸则|X7=a-T-4o c）|&经角度部分积分后其径向部分,I 1 层 Z（Z+1）”百 6Xj=9+-o 2 _ Vloc-2r dr2 2r2 r得 一 严 2口*(*1 d2 Z(Z+1)v 口(r)=/力=匕+不不/-4。丁rR 1 i(j I 1)=/drU：(r)j+而 m VLOC(r)g)这样场J-约z对应上式中第1,234项。第1项为r R（；（）出（）（与/）办=/%(?()&-Gr,dr=r 心.)显然，第3项与第4项相减之差为0 第2项为1 rR 4 1 rR J22/谩点淳2/()漉呜()分部积分后可得:=%14则Bij-岑=G-Q)6岫R ：(一)%(二)(3.16)采用同样的方法，相同结果对AE波函数0也可以得到，并且与(3.16)式相减，可以取使两 者微分在r=R时，两者微分相同，则不难得到：Bij 岑=(，ej)Qiji Qij r r(3.17)如果范数不变条件满足即。尸0,则为为厄密算符。因此根据(3.12)及(3.14)据此定义的非 局域势的算符v l也是厄密的。第三步，D.Vander bil t证明广义范数不变条件Qj=0,实际并不必要。定义广义交迭算 符S=l+Qy血%(3.18)同时，重新定义非局域势的鹰势Vnl=DQ氏为(3.19)这里。=+CjQij(3.20)则不难得到 pMr=亦血r(3.21)6iS6jr=R+3 =R+1)需 Rm,n I R=r+H Qm.n (5-+5+m,n I kR+Qmimjn-R+Qij)R R m)n不难证明R j也是厄密的Dij D%=Bij+CjQij (Bji+.Qj。*=+(与一=0故算符D也是厄密的.下而容易证明鹰波函数满足本征方程(育S)|%=015即T+VlqC+)(B/nn+nQm,n)/3m fn|iQ(1+:Qm,n/3m m,n m,nT VlOC+(3n|6 (j)i+.)2 Qm,nTTi n m,n m,n m,n上式左端显然为0,右边若用&|作用，即为：Qmnmimj&)Qm,nmi6 Tnj m,nm,n下面来证明:d 入0=pe=e.(3.22)根据前面的结果有(T+4oc+%v l)M=eS|0e 同理(T+VlqC+VL)|+A=(C+c)S|0e+Ze 上而两式分别乘以 在+/及0/，然后上两式相减，由于为厄密算符：则=一|V2+Vloc+v lIc+Ac一=6 利用格林定理：2/dm 成-e+Ae-=两边除以,取 一 0则1 f 八4*在 迎加/X iqm、LRds屹加一次加=R即；,概，*山=R Nj CuC Cy/yzg取 4=驾2可得:时姆展而也整川=河灯”(3.23a)16(3.23a)式是规范不变震势所必须满足的条件，(第4条)。因规范不变震势中。产。的限制实 际上是不必要的。这意味着每一个AE波函数。都可以独立的选择一个合适的&,而且不 必满足规范不变的条件，唯一的约束就是在r=R处所与仍需光滑连结，具有对能量的导 数相同。这样原来的一些局域性很强的价轨道，如。的2P轨道的鹰波函数也就可以选择的 尽量平缓，从而大大减少了这些鹰轨道平而波的基函数，从而对减小工作量十分有利。对 于能带为n,波矢为k的鹰波函数=nnrkk,(3.236)(3.23b)式保证了在厂“区域与实际波函数幅度相同。由于放宽了规范不变条件，在厂厂“的区域，两者并不一致，为了补偿这一部分电荷的损失，可以重新定义电荷密度 为：&(r)=T OtQnk+工 P Q/G)(3.24)n,k这里=与。疝。疝向(3.25)nkQij6=5WWW)-娘丁(3.26)这样价电子电荷密度可推出价电子总数Nv=而&=/&T味Qnk+dT Qj n,k n,k 2，J=2。/(1+ADI=WkWnk=Nvn,k J n,k因此价电子的总数保持不变$342 USPP的能量，哈密顿量体系总能量:价电子乂的总能量可以表示为区况&,必=&|一：+VNL+：/d.i 幺 2 J J r rf+Excn+I d”院n+U R (3.27)这里八为电荷密度，Exc为交换关联能。必为离子实间的相互作用，局域鹰势 中包括吟洗=立吟洗，-时以及非局域度势珠乙=汇D仅 J(3.28)qnm-j drQnmfn,m,I则基态波函数6应该是在总能量(3.27)式在=蛤的约束条件下取极值，所得到 的本征震波函数。因而有：前抵=q S仇=8&(3.29)我们下面将根据(3.27)-(3.29)来确定USP P)势下的等效哈密顿量方.从(3.27)式可以看到,其最后一项显然与电子密度无关。而第一项为6的函数，2,3,4项都是电子密度n的函数。由于n=I阳+Q：,mi不难求得:2=昕)町河+(3.30)这样总能量中的第一项：丁=+VGOi (3.31a)0cpi/总能量中的第二项：6Eh _ 8Eh 6nr 5楣 加5植 总能量中的第三项:n(rQr-rf(/)i*严 /加及Q)M(3.316)BExc 间 Hi加常胃黑二+左尸 +因 /d，Q鼠(尸)(3.31c)n,m,I这里定义从XC=Exo 间 加18为交换关联势.总能量中的第四项:Eloc _ Trion+1/不洗(，)Q用(/(3.31d)n,m,I 综合以上结果，若定义如下的屏蔽的局域有效势匕Ve)=/%+/疝卷J+4x(3.32)J r r I叫冲=口孔 十 dQ(3.33)注意,。案正将是表示P P的参数，而。上.通过匕仃依赖于波函数及电荷密度，因此在每 一个计算总能量的自洽迭代的循环中，每次都须更新。这样，有效的哈密顿量为H=-|v2+K/+E/J或 CI(3.34)n,m,I$343 USPP下动力学行为及离子力的计算在USP P下，根据CP动力学方法所得到的广义拉格朗日量，对应的应为:C+2 52 必居Etot(j)i/?/+5(3.35)这里Nj 0z,R =&杨加=。为约束条件。(3.35)中其余各量的物理含义同 前(3.2)式，则相应弓I入拉格朗日乘子之后的欧拉方程为：向=(3.36)Fi=Mi Ri=dEtot dRjdSa Ri(3.37)若体系的哈密顿量H定义如(3.34)式，离子间相互作用能量用。必表示，贝小Etot=Etot+URt(3.38)i(3.38)写成这样的形式，主要是为了把总能量中依赖于离子位置的U必一项分离出来。仔细分析比较总能量尺”的表达式，可以发现由离子运动引发离子所受力与下列几个量有 关。19，鼠=(3.39)i4m=E N (3.40)八=+QnmPnJ 乳)不难求得:詈=口舄。，+(3.42)oRi i oRj oRj*=Eazj+(3.43)a Ri zj dRj dRj9g=Q黑纯1+鸣力黑(3.44)8R!念 9R!dRi从(3.37)式可以分解得出“离子所受力”与下面三项有关J尸 V-as J 尸 dURj三=12 A牙用=薪刃 dRj dRj则Aa岫h出dRjdr-0V愣c6dRi一/d/24也受吓臀心J(yRi n,m oRj n,m oRj八一工。鼠给-E/日出务P n,m OKi n)m)I dH注意在上式导出中，（3.27）中头三项与瓦无关。而=E DnmPnm n,m,I n,m,I而DL=+/匹万如其中。黑一项仅与非局域势相关，与知无关。-dF2=a R N,。/（1+52”馆 I Bn Bm 1）1。）/n,m20d二 dRr a必m注意根据q叱的定义，它应与必无关，这样综合以上结果，离子所受力为:du dRi尸=尸1+F2+用/曲夕需-/而/E%心-E人察+的唾(3.45)n)m dRj n)m Uiii(3.45)中最后一项是约束条件的贡献，这个离子力的计算表达式虽然复杂，好在平面波基 与离子坐标无关，所以计算上还是可性的。$3.5 CP方法的具体实现图3 给出了从头计算分子动力学具体实现的基本轮廓，整个模拟过程主要通过两步实现，单 独的电子系统的动力学及电子和离子耦合系统的动力学。由t=0时刻，由密度泛函理论计 算总能量Ep,以及由Hel l mann-Feymann定理计算Ep对离子坐标瓦的微分出E,也必须知道离子初始位置下电子的单粒子轨道波函数我，所以第一步要通过电子系统的动 力学来计算原子初始构型下的基态波函数。为了给出合适的尝试波函数来进行电子虚拟动 力学模拟，可以先用少量的平面波基来展开分子轨道，按照传统的矩阵对角化技术得到分 子轨道的本征波函数，由于这时采用平面波基很小，对角化能很快完成。然后把分子轨道 本征波函数用扩展的平而波基来作为尝试波函数。有了尝试波函数，就可以按照图3的轮 廓(略去了离子的动力学过程)进行电子虚拟动力学的模拟，先计算虚拟力，哈密顿量对 每一个傅立叶函数作用然后求解波函数的运动方程由超=-Hks我，由于这一方 程没有包括正交化约束条件产生的“虚拟力所以下一步要对波函数进行正交化，得到新 的波函数。最后跳过离子运动方程的计算，只要把新的波函数作为力十力时刻的输入，就 可以进行下一步的电子动力学模拟，直至波函数收敛。第二步，得到初始本征波函数后，整个电子一离子的耦合系统就处于B。面上，就 可以开始耦合系统的动力学模拟。按照图3所示，力=0时，由第一步已知离子的初始位 置下电子的单粒子轨道基函数圾，可方便地由密度泛函理论计算总能量Ep,以及进行 电子虚拟动力学积分。通过波函数的正交化，完成电子系统的虚拟动力学过程，从而得 到力十力时刻的波函数。然后由H-F的力定理，可求得力=力十力时刻的离子坐标，21i=0个 分 子 动 力 学 循 环一对每个点-对每个电子态 I计算 M/d-工)=/C(j)FI t-一对每个平面波分量求解电子运动方程;-7二-荷正交化平面波-对每个离子|求解离子运动方程：1=_ V】E 增加一个水从头计算分子动力学计算的流程图他22倒空间I C i(g)|32门(之)Ek，g V,E 件丫Vh(G Eh 十Vf EfLVtH(g)L OC,C i(g)N-FFTFFTF F T一N-fFT_ _ _实空间 巧IP#）L（：）+vm（7）lL%不（厂）I I匕小（尸）或（厂）L 一6E/&i(g)=f t i(g图率 总能量计算的流程图，其中N.F F T表示N次F F T变换既图4中给出了模拟过程中所需要的总能Ep,HF力以及虚拟力计算的基本步骤并且标明 各自方便的计算空间。$3.5,1 NCPP下的CP方法的具体实现在NCP P下的MD模拟中总能量E可以写为：Et=Ek+Exc+Eh+Eps+U（3.46）这里其代表系统的总能量，以为电子动能。在倒空间很容易计算=|EG2C*(G)G(G)(3.47)9 J乙I小JZ/ZjRi RjU为离子之间的相互作用能。号。为交换关联能。知道了其交换关联能的密度，即可由/d走w P/）求得。Eh为电子之间的相互作用能，即Har t r ee能量。(3.48)前面已经讨论过，一旦在实空间中容易用电子密度P表示，Eps为鹰势23能与s=/dfV P e (3.49)这里的V为所有的离子作用到电子上的势能，可以表示为每个离子震势的迭加=EiVps 行一%)(3.50)在NCP P下，可以给出：Vps(r)=Vioc(r+Ez A(r)|Zm lm(3.51)(3.51)中第一项为鹰势的局域部分，第二项为鹰势的非局域部分，它是依赖于|/加态 的。这样能量Eps就可分为两部分：然(局域)卑力(非局域)=/而P e 瓦)(3.52)E尸他(3.53)应该注意的是，由于库仑力的长程性，如果单独计算E,U和或会出现不收敛的现 象，要格外小心。如果把这三项作统一处理，将来带来便利。为此对U项用电荷的密度分 布来代替点电荷，这里采用高斯(Gaus s)分布来代替点电荷分布：Pi&-总)=高户(在广一孙禺(3.54)决定了 I点的离子电荷分布密度。这样U=|Ei,jfdrf 及厂人-Eself-Eovrl(3.55)(3.55)中第一项是对所有LJ求的，并没有的限制，因此必须减掉电荷的自相互 作用项。不难求得：(3.56)心切项是由于点电荷相互作用与高斯密度分布电荷相互作用之差所引起的工 Z/Zj2 I丰j Ri-Rj(3.57)经过计算=甫为-er fc(x)为余误差函数把U,Eh，/求和得：7U+EH+Es=y 而加宣+诉Pe&VL 即一露O Eseif EoML(3.59)这里全电荷密度.=.+/,及吃c=-7山(篇)(36。)24后一项是/d/号的贡献。点电荷用高斯密度分布代替后，就可以象电子的电荷分布一样在倒空间中计算。从 而在倒空间中这三项都可以统一处理。把电子电荷与离