密度泛函理论与从头计算分子动力学
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
密度泛函理论与从头计算分子动力学
$ 1 引言
此处有误 自从上世纪60年代以来,密度泛函理论(DFT)建立.并在局域密度近似(LDA)下导出著名
的Kohn-Sham(KS)方程以来,DFT一直是凝聚态物理领域计算电子结构及其特性的最有力 的工具. 近几年来DFT与分子动力学相结合,在材料设计,合成,模拟计算和评价诸多方面有 明显进展,成为计算材料科学的重要基础和核心技术. 特别在量子化学计算领域,1987年以 前主要用Hartree-Fock(HF)方法.但近年来,用DFT的工作以指数增加.以致于HF方法应用已 经相当少.W.Kohn因提出DFT获得1998年诺贝尔化学奖. 已经表明了DFT在计算化学领域 的核心作用与应用的广泛性. 分子动力学计算机模拟是研究复杂的凝聚态系统的有力工具.这一技术既能得到原 子的运动轨迹,还能像做实验一样作各种观察.对于平衡系统,可以在一个分子动力学观 察(observation time)内作时间平均来计算一个物理量的统计平均值.对一个非平衡系统过 程,只要发生在一个分子动力学观察时间内(1-10ps)的物理现象也可以用分子动力学计算进 行直接模拟.可见数值实验对理论与实验的有力补充,特别是许多与原子有关的微观细节,在 实验中无法获得而在计算机模拟中可以方便的得到. DFT与分子动力学(MD)相结合可以有大量不同类型的应用. 如晶格生长,外延生长,离 子移植,缺陷运动,无序结构,表面与界面的重构,电离势的计算,振动谱的研究,化学反应的 问题,生物分子的结构,催化活性位置的特性以及材料电子结构和几何结构.固,液体的相变 等.现在这些方法已发展成为成熟的计算方法.DFT 的另一个特点是,它提供了第一性原 理(first-principal)或称为从头计算的理论框架.在这个框架下可以发展各种各样的能带计 算方法.虽然在DFT的所有这些实际应用中,几乎都采用局域密度近似(LDA),这是一种不能 控制精度的近似,因而DFT方法的有效性在很大程度上要看结果与实验的一致性.人们没 有任何直接的方法可以改善LDA的精度. 然而DFT允许发展别的方法加以补充.例如广义 梯度近似(GGA)等方法,把密度分布n(r )的空间变化包括在方法之中,实现了比较大幅度减 少LDA误差的目的. 相比较传统的量子化学方法,如组态相互作用(CI)方法,DFT+MD方法显然可以用于数 百个原子的大分子体系.但对于具有强关联相互作用的体系,似乎化学家更愿意应用CI方 法.而对于凝聚态物理领域,DFT+MD方法可以相当精确的计算材料的电子结构和相应的 许多物性. 在DFT获得巨大成功的同时,也有些不可忽视的弱点与困难.针对这些问题已经发展了许
V
T,V为严格的动能及电子间的势能函数。T0 为无相互作用的电子气动能，VH 为电子间 的库仑能，VX 为电子间交换能。因而，不难得到上述两种体系下能量差别仅仅在于电子 相关能VC . 则 VC =Ee − EHF =T-T0 则HK的普适函数FHK 可以写为： FHK = T +V +T0 −T0 = T0 +V +(T − T0 ) = T0 +V +VC +VH −VH = T0 +VH +VC +(V − VH )
6 波函数φi 上只需O(NM)次操作。势能Vef f 是局域的，则作用到波函数φi 上也可方便的用 快速付立叶变换技术(FFT)在实空间计算, 也只要O(NMln M )次操作。因而工作量以线性 或O(NMln M )次增长,这与传统的标准对角化技术相比,工作量有显著减小.
$ 2.4 局域密度近似下的离子间相互作用势的计算
xc (r )只取决于该点的ρ(r )
,整个系统的交换关联能为
xc [ρ(r )]
EXC =
drρ(r)
(2.13)
至 于 交 换 关 联 势EXC 的 具 体 表 达 式 不 断 有 著 者 给 出 。 目 前 用 的 最 广 泛 的 是Cepeley和Alder在 八 十 年 代 用Monte Carlo方 法 导 出 的 解 析 表 示(Phys.Rev.lett(566568)Vol.45.No.7)。从方法推导近似来看，LDA只局限于那些电子密度缓变的系统。因 而，在1996年，Perdew,Burke,Ernzerhof把它推广到包括电子梯度密度在内的函数，并给 出了解析表示式。称为GGA方法。（Generalized Gradient Approximation.Phys.Rev.lett Vol.77.No.8(3895-3897)),如包含电子自旋极化(LSD)情况在内，则（2.13)可写为
$ 2 基本理论
首先，对包括诸多电子，离子，原子的凝聚态体系来说，面临的首要问题是，由于原 子核较重而且运动中无法跟随电子的运动。因而可以把电子与原子分开考虑，原子实作为
3 带有正电荷的外场存在。这就是所谓的绝热近似或称BO近似（Born–Oppenheimer).因而 我们可以把电子的哈密顿量写成 ˆ =T ˆ+V ˆ +V ˆext H (2.1)
(2.9)–(2.12)式就构成了KS方法的基本方程组。 遗憾的是，它还是形式上的东西，因为包含电子多体效应交换关联效应的EXC 的具体 表达式还不清楚，需作某种近似.现在广泛应用的是局域密度近似（LDA），即把非均匀 电子系统分割成一些小块。在这些小块中认为电子分布是均匀的。这样r处的子块中的交 换关联能密度
i
(2.16)
δE 上式中的点· 代表对时间的微分。变分 δφ ∗ = HKS φi ,依赖于φi ，求解上式时，为保证波
函数的正交性，还得加上一个约束条件，虚拟动力学方程是： ˙ i (r, t) = − φ δE δφ∗ i (r, t) +
j
∧ij φj
(2.17)
实际上，先用无约束的动力学方程(2.16)计算φi ，然后通过适当的正交化过程对φi 进行 ˙ i = 0 ,即HKS φi = j ∧ij φj ,把相应的约束 修正，并得到系数∧ij 。当系统达到极小值时，φ 矩∧ij 对角化后，既可以得到KS方程的本征矢与本征值。经验表明，对于一定的体系坐 标E只有一个极小值。这样速降法也能得到能量极小值。这种方法与传统的方法相比省去 了大量的对角化计算。 速降法的有效性取决于波函数达到收敛的有效步数，这个步数可能很多，特别是在低 对称下，人们也做了许多改善速降法的尝试，如有更复杂的积分方程，以及改善的速降 法，及共扼梯度法等等。 ˆ KS 作用到φi 的每步计算中都需要适当 象速降法及共扼梯度法这样的算法在哈密顿量H 的正交化操作。如果φi 是用平面波展开的（设N个单电子轨道均用M个平面波展开）， 则解HKS φi 就需要NM 2 次浮点计算。加上正交化需要N 2 M次操作,计算量很大。为了减少 计算量，可以把算符HKS 分成动能项与势能项，动能项在倒格矢空间是对角的，作用到
$ 2.2 Kohn-Sham方程
K―S方程诞生于1965年,K―S方程提出后，才真正使DFT成为实际计算可能。设总能 量函数Ee [ρ] (包括电子间相关能）及在Hartree-Fock近似下的总能量函数为EHF [ρ]（不包 括电子间相关能） 则 Ee =T+V (2.5) EHF = T0 + (VH + VX ) (2.6)
$ 2.3 Kohn-Sham方程的求解
在传统的固体，分子的局域密度泛函电子结构计算中，KS自洽方程组一般可以这样求 解：先给定一个电荷ρo (r),比如说可以选用各原子电荷密度的迭加，然后代入(2.12)计算有 ˆ KS 的本征矢。即φi (r)。然后可以求得新的分 效势Vef f ,再通过对角化来求解(2.10) ，得到H 布电荷密度ρo (r)，继续进行迭代直致自洽。由于对角化计算工作量以M 3 增加(M为KS轨 道展开所需基函数的个数)，因而对有数百个原子的系统工作量非常庞大。 KS方程组也可以通过其它方法解。为了求得能量泛函对电子密度变分的最小值，我们 可以在电子自由度φi 的空间中引入一个虚拟的动力学过程，来使泛函达到极小值。最简单 动力学考虑是速降法。 ˙ i (r, t) = − φ δE δφ∗ i (r, t)
VC VX
(2.7)
4 =T0 +VH +(VX + VC )
VXC
(2.8)
这里VXC 被称为交换–关联势函数，这样总能量泛函可以写为 EVext [ρ] = T0 [ρ0 ] + (1/2) 这里交换关联势可由下式给出 ˆXC = δEXC [ρ] V δρ 对应的系统哈密顿量（通常称为KS哈密顿量）写为 ˆ ks φi = [− 1 ∇2 + V ˆef f ]φi = Ei φi H 2 i 这里单粒子波函数满足 ρ( r ) =
$ 2.1 Hohenberg-Kohn定理
1964年，HK提出HK两个定理： ˆ的 1) 多电子系统在外场势Vext 作用下，其基态的电子密度ρ(r), 与基态下任意力学量O 可观测量是一一对应的，且是严格的基态电子密度ρ(r)的泛函. ˆ |Ψ >=O[ρ] < Ψ|O (2.2) ˆ 为哈密顿量H ˆ ，则基态的总能量函数H[ρ]=EV ext [ρ]具有如下形式 2) 若O ˆ+V ˆ |Ψ >+< Ψ|V ˆext |Ψ > EV ext [ρ]=< Ψ|T (2.3) 而对任意多电子体系的普适量FHK [ρ]为 ˆ+V ˆ |Ψ > FHK [ρ]=< Ψ|T EV ext [ρ]=FHK [ρ]+ ρ(r)Vext (r)dr 这里，我们不准备证明上述两个定理，只说明如下3点 1）任何ρ与力学量或势之间存在一一对应性（one to one correspondence) 2)普适函数FHK [ρ]可以很容易用密度算符表示。 3)HK第二定理可以使我们利用变分定理，通过求能量极小值来获得基态的电子密度。 (2.4)
传统的,如果给定原子坐标{ RI },由从头计算方法推导离子间相互作用势V[{ RI }],然后 进行第一性原理的分子动力学模拟,一般有三个基本步骤,即在每个MD 过程中:(1)对于给定 的{ RI }求解自洽的KS方程.(2)根据Hellman–Feymam的力定理计算每个离子受到的力.(3) 求解牛顿运动方程 ¨ I = −∇R V MI R I 但需要同时对电子和离子的自由度进行优化(Bendt,Zunger).利用速降法,通过下面的方 程组就可以简单的实现这样的想法. ˙ i = − δE + µi φ δφ∗ i ∧ij φj
i=1 N
drdr
ρo (r)ρo (r ) + |r − r |
Vext [ρ]dr + EXC [ρo ] (2.9)
(2.10)
φ∗ i (r )φi (r ) dr
(2.11) (2.12)
ˆef f = Vext (r) + V
ρ(r ) δEXC [ρ] + δρ(r) |r − r |
2
多不同的方法,这些方法可以用Kohn-Sham方程的有效Hamiltonian的各部分和波函数构造 上的考虑进行归类,如图一所示 传统的DFT―LDA是预言多电子体系基态性质的理论。对于激发态性质的描述总是与 实验不符合，这固然因为（1）激发态本身存在着复杂性质。（2）DFT―LDA理论存在 着对激发态描述困难。此外，对于具有强关联相互作用体系LDA理论描述不够好。因而 发展出LDA+U和LDA++等方法。此外关于大原子数的复杂体系，近几年来发展了各种 线性标度的方法，也称为O（N）算法。为研究复杂体系提供了有力的工具。
LSD EXC [n ↑, n ↓] =
d3 r n
unif xc [n
百度文库
↑, n ↓]
(2.14)
而在GGA下，则为
LSD [n ↑, n ↓] = EXC
d3 r f (n ↑, n ↓, ∇n ↑, ∇n ↓)
(2.15)
5 现在GGA方法已经得到广泛的应用。总之，局域密度泛函理论，自上世纪60年代提出以 来，经过许多人的发展与完善，已广泛的应用于各领域多原子体系的电子结构计算。由于 在此理论近似下，单电子运动方程中的交换关联势形式简单，使计算工作量大大减少。因 此，现在大多数能带或者原子集团模型中的计算方法，都是在局域密度泛函理论近似的基 础上建立起来的。
$ 1 引言
此处有误 自从上世纪60年代以来,密度泛函理论(DFT)建立.并在局域密度近似(LDA)下导出著名
的Kohn-Sham(KS)方程以来,DFT一直是凝聚态物理领域计算电子结构及其特性的最有力 的工具. 近几年来DFT与分子动力学相结合,在材料设计,合成,模拟计算和评价诸多方面有 明显进展,成为计算材料科学的重要基础和核心技术. 特别在量子化学计算领域,1987年以 前主要用Hartree-Fock(HF)方法.但近年来,用DFT的工作以指数增加.以致于HF方法应用已 经相当少.W.Kohn因提出DFT获得1998年诺贝尔化学奖. 已经表明了DFT在计算化学领域 的核心作用与应用的广泛性. 分子动力学计算机模拟是研究复杂的凝聚态系统的有力工具.这一技术既能得到原 子的运动轨迹,还能像做实验一样作各种观察.对于平衡系统,可以在一个分子动力学观 察(observation time)内作时间平均来计算一个物理量的统计平均值.对一个非平衡系统过 程,只要发生在一个分子动力学观察时间内(1-10ps)的物理现象也可以用分子动力学计算进 行直接模拟.可见数值实验对理论与实验的有力补充,特别是许多与原子有关的微观细节,在 实验中无法获得而在计算机模拟中可以方便的得到. DFT与分子动力学(MD)相结合可以有大量不同类型的应用. 如晶格生长,外延生长,离 子移植,缺陷运动,无序结构,表面与界面的重构,电离势的计算,振动谱的研究,化学反应的 问题,生物分子的结构,催化活性位置的特性以及材料电子结构和几何结构.固,液体的相变 等.现在这些方法已发展成为成熟的计算方法.DFT 的另一个特点是,它提供了第一性原 理(first-principal)或称为从头计算的理论框架.在这个框架下可以发展各种各样的能带计 算方法.虽然在DFT的所有这些实际应用中,几乎都采用局域密度近似(LDA),这是一种不能 控制精度的近似,因而DFT方法的有效性在很大程度上要看结果与实验的一致性.人们没 有任何直接的方法可以改善LDA的精度. 然而DFT允许发展别的方法加以补充.例如广义 梯度近似(GGA)等方法,把密度分布n(r )的空间变化包括在方法之中,实现了比较大幅度减 少LDA误差的目的. 相比较传统的量子化学方法,如组态相互作用(CI)方法,DFT+MD方法显然可以用于数 百个原子的大分子体系.但对于具有强关联相互作用的体系,似乎化学家更愿意应用CI方 法.而对于凝聚态物理领域,DFT+MD方法可以相当精确的计算材料的电子结构和相应的 许多物性. 在DFT获得巨大成功的同时,也有些不可忽视的弱点与困难.针对这些问题已经发展了许
V
T,V为严格的动能及电子间的势能函数。T0 为无相互作用的电子气动能，VH 为电子间 的库仑能，VX 为电子间交换能。因而，不难得到上述两种体系下能量差别仅仅在于电子 相关能VC . 则 VC =Ee − EHF =T-T0 则HK的普适函数FHK 可以写为： FHK = T +V +T0 −T0 = T0 +V +(T − T0 ) = T0 +V +VC +VH −VH = T0 +VH +VC +(V − VH )
6 波函数φi 上只需O(NM)次操作。势能Vef f 是局域的，则作用到波函数φi 上也可方便的用 快速付立叶变换技术(FFT)在实空间计算, 也只要O(NMln M )次操作。因而工作量以线性 或O(NMln M )次增长,这与传统的标准对角化技术相比,工作量有显著减小.
$ 2.4 局域密度近似下的离子间相互作用势的计算
xc (r )只取决于该点的ρ(r )
,整个系统的交换关联能为
xc [ρ(r )]
EXC =
drρ(r)
(2.13)
至 于 交 换 关 联 势EXC 的 具 体 表 达 式 不 断 有 著 者 给 出 。 目 前 用 的 最 广 泛 的 是Cepeley和Alder在 八 十 年 代 用Monte Carlo方 法 导 出 的 解 析 表 示(Phys.Rev.lett(566568)Vol.45.No.7)。从方法推导近似来看，LDA只局限于那些电子密度缓变的系统。因 而，在1996年，Perdew,Burke,Ernzerhof把它推广到包括电子梯度密度在内的函数，并给 出了解析表示式。称为GGA方法。（Generalized Gradient Approximation.Phys.Rev.lett Vol.77.No.8(3895-3897)),如包含电子自旋极化(LSD)情况在内，则（2.13)可写为
$ 2 基本理论
首先，对包括诸多电子，离子，原子的凝聚态体系来说，面临的首要问题是，由于原 子核较重而且运动中无法跟随电子的运动。因而可以把电子与原子分开考虑，原子实作为
3 带有正电荷的外场存在。这就是所谓的绝热近似或称BO近似（Born–Oppenheimer).因而 我们可以把电子的哈密顿量写成 ˆ =T ˆ+V ˆ +V ˆext H (2.1)
(2.9)–(2.12)式就构成了KS方法的基本方程组。 遗憾的是，它还是形式上的东西，因为包含电子多体效应交换关联效应的EXC 的具体 表达式还不清楚，需作某种近似.现在广泛应用的是局域密度近似（LDA），即把非均匀 电子系统分割成一些小块。在这些小块中认为电子分布是均匀的。这样r处的子块中的交 换关联能密度
i
(2.16)
δE 上式中的点· 代表对时间的微分。变分 δφ ∗ = HKS φi ,依赖于φi ，求解上式时，为保证波
函数的正交性，还得加上一个约束条件，虚拟动力学方程是： ˙ i (r, t) = − φ δE δφ∗ i (r, t) +
j
∧ij φj
(2.17)
实际上，先用无约束的动力学方程(2.16)计算φi ，然后通过适当的正交化过程对φi 进行 ˙ i = 0 ,即HKS φi = j ∧ij φj ,把相应的约束 修正，并得到系数∧ij 。当系统达到极小值时，φ 矩∧ij 对角化后，既可以得到KS方程的本征矢与本征值。经验表明，对于一定的体系坐 标E只有一个极小值。这样速降法也能得到能量极小值。这种方法与传统的方法相比省去 了大量的对角化计算。 速降法的有效性取决于波函数达到收敛的有效步数，这个步数可能很多，特别是在低 对称下，人们也做了许多改善速降法的尝试，如有更复杂的积分方程，以及改善的速降 法，及共扼梯度法等等。 ˆ KS 作用到φi 的每步计算中都需要适当 象速降法及共扼梯度法这样的算法在哈密顿量H 的正交化操作。如果φi 是用平面波展开的（设N个单电子轨道均用M个平面波展开）， 则解HKS φi 就需要NM 2 次浮点计算。加上正交化需要N 2 M次操作,计算量很大。为了减少 计算量，可以把算符HKS 分成动能项与势能项，动能项在倒格矢空间是对角的，作用到
$ 2.2 Kohn-Sham方程
K―S方程诞生于1965年,K―S方程提出后，才真正使DFT成为实际计算可能。设总能 量函数Ee [ρ] (包括电子间相关能）及在Hartree-Fock近似下的总能量函数为EHF [ρ]（不包 括电子间相关能） 则 Ee =T+V (2.5) EHF = T0 + (VH + VX ) (2.6)
$ 2.3 Kohn-Sham方程的求解
在传统的固体，分子的局域密度泛函电子结构计算中，KS自洽方程组一般可以这样求 解：先给定一个电荷ρo (r),比如说可以选用各原子电荷密度的迭加，然后代入(2.12)计算有 ˆ KS 的本征矢。即φi (r)。然后可以求得新的分 效势Vef f ,再通过对角化来求解(2.10) ，得到H 布电荷密度ρo (r)，继续进行迭代直致自洽。由于对角化计算工作量以M 3 增加(M为KS轨 道展开所需基函数的个数)，因而对有数百个原子的系统工作量非常庞大。 KS方程组也可以通过其它方法解。为了求得能量泛函对电子密度变分的最小值，我们 可以在电子自由度φi 的空间中引入一个虚拟的动力学过程，来使泛函达到极小值。最简单 动力学考虑是速降法。 ˙ i (r, t) = − φ δE δφ∗ i (r, t)
VC VX
(2.7)
4 =T0 +VH +(VX + VC )
VXC
(2.8)
这里VXC 被称为交换–关联势函数，这样总能量泛函可以写为 EVext [ρ] = T0 [ρ0 ] + (1/2) 这里交换关联势可由下式给出 ˆXC = δEXC [ρ] V δρ 对应的系统哈密顿量（通常称为KS哈密顿量）写为 ˆ ks φi = [− 1 ∇2 + V ˆef f ]φi = Ei φi H 2 i 这里单粒子波函数满足 ρ( r ) =
$ 2.1 Hohenberg-Kohn定理
1964年，HK提出HK两个定理： ˆ的 1) 多电子系统在外场势Vext 作用下，其基态的电子密度ρ(r), 与基态下任意力学量O 可观测量是一一对应的，且是严格的基态电子密度ρ(r)的泛函. ˆ |Ψ >=O[ρ] < Ψ|O (2.2) ˆ 为哈密顿量H ˆ ，则基态的总能量函数H[ρ]=EV ext [ρ]具有如下形式 2) 若O ˆ+V ˆ |Ψ >+< Ψ|V ˆext |Ψ > EV ext [ρ]=< Ψ|T (2.3) 而对任意多电子体系的普适量FHK [ρ]为 ˆ+V ˆ |Ψ > FHK [ρ]=< Ψ|T EV ext [ρ]=FHK [ρ]+ ρ(r)Vext (r)dr 这里，我们不准备证明上述两个定理，只说明如下3点 1）任何ρ与力学量或势之间存在一一对应性（one to one correspondence) 2)普适函数FHK [ρ]可以很容易用密度算符表示。 3)HK第二定理可以使我们利用变分定理，通过求能量极小值来获得基态的电子密度。 (2.4)
传统的,如果给定原子坐标{ RI },由从头计算方法推导离子间相互作用势V[{ RI }],然后 进行第一性原理的分子动力学模拟,一般有三个基本步骤,即在每个MD 过程中:(1)对于给定 的{ RI }求解自洽的KS方程.(2)根据Hellman–Feymam的力定理计算每个离子受到的力.(3) 求解牛顿运动方程 ¨ I = −∇R V MI R I 但需要同时对电子和离子的自由度进行优化(Bendt,Zunger).利用速降法,通过下面的方 程组就可以简单的实现这样的想法. ˙ i = − δE + µi φ δφ∗ i ∧ij φj
i=1 N
drdr
ρo (r)ρo (r ) + |r − r |
Vext [ρ]dr + EXC [ρo ] (2.9)
(2.10)
φ∗ i (r )φi (r ) dr
(2.11) (2.12)
ˆef f = Vext (r) + V
ρ(r ) δEXC [ρ] + δρ(r) |r − r |
2
多不同的方法,这些方法可以用Kohn-Sham方程的有效Hamiltonian的各部分和波函数构造 上的考虑进行归类,如图一所示 传统的DFT―LDA是预言多电子体系基态性质的理论。对于激发态性质的描述总是与 实验不符合，这固然因为（1）激发态本身存在着复杂性质。（2）DFT―LDA理论存在 着对激发态描述困难。此外，对于具有强关联相互作用体系LDA理论描述不够好。因而 发展出LDA+U和LDA++等方法。此外关于大原子数的复杂体系，近几年来发展了各种 线性标度的方法，也称为O（N）算法。为研究复杂体系提供了有力的工具。
LSD EXC [n ↑, n ↓] =
d3 r n
unif xc [n
百度文库
↑, n ↓]
(2.14)
而在GGA下，则为
LSD [n ↑, n ↓] = EXC
d3 r f (n ↑, n ↓, ∇n ↑, ∇n ↓)
(2.15)
5 现在GGA方法已经得到广泛的应用。总之，局域密度泛函理论，自上世纪60年代提出以 来，经过许多人的发展与完善，已广泛的应用于各领域多原子体系的电子结构计算。由于 在此理论近似下，单电子运动方程中的交换关联势形式简单，使计算工作量大大减少。因 此，现在大多数能带或者原子集团模型中的计算方法，都是在局域密度泛函理论近似的基 础上建立起来的。