密度泛函理论

密度泛函理论及基本应用
凝聚态物理陈阿海2010210602
摘要：本文讨论了密度泛函理论的基本原理和方法，并讨论了玻色子系统中处于谐振子势以及高斯形势垒下的两粒子问题。

随着成对的高斯势垒分别向两边移动，粒子数分布由原来被高斯势垒劈裂的两个峰过度到单个峰。

关键词：密度泛函理论，玻色子，谐振势，高斯势垒
Abstract: The basic method of density functional theory is discussed as well as a case of two bosons trapped in a harmonic external potential with Gaussian potential barrier. With the move of a pair of Gaussian potential barrier to the two edges respectively, the distribution of atoms changed gradually from the two peaks divided by the Gaussian potential barrier in the centre of the system at first to single peak again.
Key words: density functional theory, boson, harmonic external potential, Gaussian potential barrier
一、引言
密度泛函理论（density functional theory，DFT）是当前研究量子多体系统的重要理论之一，在物理与化学领域均有广泛的应用。

DFT的基本出发点是将量子多体系统的性质全部表达为密度的泛函，从而达到简化系统的目的。

密度泛函理论的基本概念最早起源于十九世纪二十年代的Thomas-Fermi模型（TF）[1][2]，但其真正意义上的理论开始于六十年代Hohenberg与Kohn的研究，其在1964年发表的著名论文[3]奠定了密度泛函理论的基础。

其后，Kohn与Sham在1965年发表的论文将DFT理论推向了实际应用的水平[4]，其提出的Kohn-Sham方法成为密度泛函理论的基本应用框架。

Kohn凭借其对密度泛函理论的贡献获得了1998年的诺贝尔化学奖。

在其后的几十年间，DFT在多个领域得到了进一步的发展。

在量子多体理论的研究中，DFT已发展成为重要的数值计算手段之一，如同量子蒙特卡罗（quantum monte carlo, QMC），密度重整化群（density matrix renormalization group, DMRG），精确对角化（exact diagonalization, ED）等数值方法一样，是当前量子多体系统最主流的数值计算方法之一。

为了进一步扩大DFT的应用范围，相继出现了处理含时问题的含时密度泛函理论（time dependent density functional theory, TDDFT）[5-10]，处理自旋系统的自旋密度泛函理论（spin density functional theory, SDFT）[11-17]，处理有限温度问题的有限温度密度泛函理论（finite temperature density functional theory, FTDFT）[18-20]，整合含时密度泛函理论和自旋密度泛函理论的含时自旋密度泛函理论（time-dependent spin-density functional theory, TD-SDFT）[21]等等各种衍生理论和方法。

本文讨论DFT理论的基本原理和方法，进一步利用一个简单的例子讨论DFT的实际应用。

在该例子中，主要讨论受谐振势和高斯型势垒作用的一维玻色系统两粒子问题，其中高斯势垒由势阱中央出现并成对向系统两侧移动。

主要讨论了该外势作用下系统粒子分布的变化特点。

本文组织如下：第二节中给出了DFT理论的基本框架和特点，主要包括Hohenberg-Kohn 定理，Kohn-Sham方程以及密度泛函理论计算的主要特点；第三节中讨论了零温下受外势作用的一般一维玻色系统的DFT；第四节具体给出了上述所提的两粒子问题的计算结果并作了简单的讨论；最后给出简单的小结。

二、DFT理论的基本框架和特点
基本的DFT理论基于两大基本原理，即Hohenberg-Kohn定理。

定理1指出量子系统的外势与系统的密度分布一一对应（可相差一个常数），或者说，外势是系统密度的泛函，而这种映射是双向的，即系统外势与系统密度分布具有相互确定的一一对应关系。

定理2即系统密度的变分原理，其指出给定有关粒子数和外势的系统密度分布，那么，使得系统能量最小的密度分布即是系统的基态密度分布，相应的能量即为系统的基态能量，即：
int 0[()][()][()][()]ext E n r T n r U n r U n r E =++≥ (1)
仅当0()()n r n r =时取得等号，其中[()]T n r 为系统动能，[()]ext U n r 为外势所对应的能量，int [()]U n r 为粒子相互作用的能量，()n r 为任意密度分布，0()n r 为基态密度分布，0E 为基态能量。

基于Hohenberg-Kohn 定理，Kohn-Sham 等式提供了实际应用的框架。

以下为Kohn-Sham 等式的简单推导。

从密度泛函的角度出发，量子系统的能量可表为
[()][()][()][()][()]ext ncl E n r T n r C n r U n r E n r =+++ (2)
其中，[()]T n r 为系统动能，[()]ext U n r 为外势相应的能量，[()]ncl E n r 表示所有非经典部分的能量，[()]C n r 为Hartree 能（库仑相互作用能量）：
2()(')[()]'2|'|
q n r n r C n r drdr r r =-⎰ (3) 将动能表示如下两部分：
[()][()][()]S C T n r T n r T n r =+
其中[()]S T n r 为类似Hartree-Fock 处理的单粒子近似：
21[()]||2S i i i
T n r ψψ=-<∇>∑ [()]C T n r 为额外项。

带入原式并将[()]C T n r 与[()]ncl E n r 合并为[()]XC E n r ，统称为交换关联能，得：
[()][()][()][()][()]S ext XC E n r T n r C n r V n r E n r =+++ (4)
根据原理二，利用
()()0n r dr N n r dr δ==⎰⎰，
引入拉普拉斯乘子μ，取变分： [[()]()]0()
E n r n r dr n r δμδ-=⎰ (5) 可得： [()]U [()][()][()]()()()()
S ext XC T n r n r E n r C n r n r n r n r n r δδδδμδδδδ+++= (6) 令 H [()]()()[()]ks xc V n r r r n r υυυ=++ (7)
称为Kohn-Sham 势（KS 势），其中 U [()][()][()]()()[()]()()()
ext XC H xc n r E n r C n r r r n r n r n r n r δδδυυυδδδ===， ， 分别为Hartree 势，系统所受外势，交换关联势。

则原式为：
[()][()]()S ks T n r V n r n r δμδ+= 最终得到Kohn-Sham 方程：
22[()]()()2ks i i i V n r r r m ψεψ⎛⎫-∇+= ⎪⎝⎭
(8) 得到Kohn-Sham 方程，即可利用自洽场近似来求解系统基态。

根据以上的推导可见，基本的DFT 理论主要处理系统基态的问题，尤其是基态密度分布和基态能量。

对于求解Kohn-Sham 方程，H ()r υ与()r υ均易得，求解的难点在于交换关联势[()]xc n r υ，通常采用局域密度近似（local density approximation ，LDA ）来处理，随后也有一些问题借用其他数值方法，比如QMC 、DMRG 等。

一种常用的处理方式为：
hom 0[()]c H x n r υμυμ=--
其中hom μ为均匀系统化学势，H υ为Hartree 势，0μ为无相互作用系统的化学势。

采用密度泛函理论计算时存在一定的局限性，比如通常无法得到系统波函数，无法计算系统关联性质等，同时其计算的精度也无法达到非常高，其适用的系统也有一定的限制（比如无法很好的处理具有陡变外势的问题）。

然而密度泛函理论有其自身优势和特点，作为典型的第一性原理计算，对于其适用的系统，具有极好的可靠性和运行效率，利用密度泛函理论可以计算上千个粒子的系统并保持良好的精度，其计算速度也极为优秀。

三、一维玻色系统DFT （T=0）
如上所述，求解Kohn-Sham 方程主要难点在于构造交换关联势[()]xc n r υ，采用LDA 处理时，利用均匀系统的化学势代替非均匀系统的化学势：
hom hom '0[]=[']|xc n H n n n υυμυμ→=-- (9)
对于玻色系统，无相互作用系统的化学势0μ=0，而H υ易得。

问题归结于求解均匀系统化学势hom μ。

hom μ的求解可依据均匀系统的基态精确解。

一维δ排斥相互作用的均匀玻色系统的精确解最早由[22] Lieb 和Liniger 给出，采用了所谓的Lieb-Liniger 模型：
2
2
1D 2+g ()2i i i i j
i H x x m x δ>∂=--∂∑∑ (10) 分别为动能量和相互作用项，其中1D g 为排斥相互作用强度。

在周期边界条件下，可利用Bethe-ansatz 得到其精确解，其基态能量密度可表为：
22()()2n n e m εγ=
(11) 其中1122()D D mg mg g nh n h
g γ==
=，为无量纲量，n 为系统密度。

在热力学极限下，()e γ满足 321
))((Q Q e k n k dk γρ+-=⎰ (12)
又()k ρ与Q 分别满足： 22112()()'22(')Q
Q
g k k dk g k k ρρππ+--++=⎰ (13) ()Q
Q k dk n ρ+-=⎰ (14)
令k QX =，求解该耦合方程，易得()e γ。

根据文献[23]，()e γ可通过下式得到很精确的拟合：
221223122/(1)3()1p p e q q p γγγπγγγγ
+=++++ (15)
其中2121=2.05737, =0.0688097, =3.16861, =0.898297p p q q ，本文采用该拟合方程。

则均匀系统的化学势为
22([])2()2hom e h n e m γμ
γγγ∂-∂= (16) 需求解()e γγ
∂∂，对(11)求导并求解耦合方程，或直接对(15)式求导可得，此处从略。

至此，已求得所有相应项的表达式，交换关联势可表为：
hom [()]c H x n r υμυ=- (17)
带入Kohn-Sham 势(7)式，可得
hom hom H []k H s V n υυμυυμ=++-=+ (18)
带入Kohn-Sham 方程得：
22hom ()()2i i i r r m υμψεψ⎛⎫-∇++= ⎪⎝⎭
(19) 求解方程可得系统基态密度分布及基态能量。

四、具有谐振势以及高斯势垒的两粒子玻色系统
在(19)式中引入外势表达式，可求得特定外势下的系统特点，相应的密度分布为：
20()|()|n r N r ψ= (20)
此处我们考虑处于谐振势和高斯型势垒下的两粒子问题，所加的外势参看图1。

外势由一个谐振势阱和高斯势垒组成，高斯型势垒可分为两个并对称地向两边移动，可表示如下：
222()()1()()
2x s x s x x a e e υ-+--=++ (21)
其中，谐振势中已令1ω=，a 与s 为高斯势垒的参数，前者控制势垒高度，后者控制势垒
位置，当存在两个势垒时（如图1(c)(d)），0s >，令10a =，当存在单个势垒时，0s =，令5a =，从而保证外势如图1所示。

图 1. 所加系统外势的概念图。

图(a)为简单谐振势，图(b)为高斯势垒居于
系统中央，图(c)(d)为向两边移动的高斯势垒。

根据上述方法将外势带入Kohn-Sham 方程：
2222hom 2()()1()()()()22x i i i s x s i x x a e e x x m μψψεψ-+--⎛⎫⎛⎫-∇++++= ⎪ ⎪⎝⎭⎝⎭
(21) 求解方程可得系统密度分布。

控制的参数主要有粒子间相互作用强度1D g ，高斯势垒的位置和高度，系统长度的选取总是保证所有粒子束缚于势阱中并在边缘平滑趋于零。

由于
图 2. 单纯谐振外势下的两粒子分布，粒子集中于系统中央。

该系统有精确解，可用于可靠性验证[5]图 3. 高斯势垒处于势阱中心时粒子分布，势垒的存在使粒子对称分布于势垒两边。

图4. s=1时粒子密度分布。

图中红色虚线标出了两个势垒的峰所在的位置（下同）。

粒子的分布也起于该位置的边缘并对称分布于两边图5. s=1.25时粒子密度分布。

红色虚线标出了两个势垒的峰所在的位置。

由于粒子向两边移动所导致的能量增加使得粒子越过势垒重新集中于势阱中央。

从图可见红线所标识的势垒峰的位置上粒子密度不为零。

图 6. s=2时粒子密度分布。

随着势垒继续向两边移动，最终又回到如图1中的单个峰的情况，由于受到两侧势垒的挤压，其峰值要高于图1中的情况。

图7. s=3时粒子密度分布。

随着势垒的继续移动，其对处于系统中心的粒子的影响逐渐减小，粒子分布峰值又逐渐减小
图8. 高斯势垒的变化。

DFT 对于处理陡变的外势作用存在一定的局限，对于变化剧烈的外势，常常不能得到较好的结果，因此，此处不讨论高斯势垒高度的变化（总是使势垒的高度控制在 ），主要从粒子高斯势垒的位置方面讨论。

由于DFT 比较善于处理相互作用比较小的系统[24]，令0.1g =，相应计算参数设置如下：0a =；5a =，0s =；10a =，1,1.25,2,3s =，相应结果见图2-9。

图2为单纯谐振势下的两粒子分布0a =。

此时系统有精确解，可参考文献[24]，这也是对该方法有效性的一个验证，此处从略。

系统的密度分布体现了谐振势的束缚作用，同时由于较弱的粒子间相互作用，粒子集中分布于势阱中心。

图3为高斯势垒处于势阱中心的情况，5a =，0s =。

从图可见，势垒的存在使粒子对称分布于势垒两边。

图4中，10a =，1s =，势垒分裂为两个，并处于1i =±的位置。

此时势垒在0i =交错，图中红色虚线标出了两个势垒的峰所在的位置。

粒子的分布也起于该位置的边缘并对称分布于两边，类似于图 3 的情况，但分离得更远。

从图2至图4，由于处于系统中央区域的势垒和谐振外势的作用使得粒子分布分裂为两个分并随着势垒的移动被向两边挤压，此时系统的能量不断增加，当粒子分布移动所引起的谐振势相应能量的增加大于高斯势垒时，将出现图5所示的临界情况。

由于粒子向两边移动所导致的能量增加使得粒子越过势垒重新集中于势阱中央，甚至有一部分几率分布于势垒峰所在位置，从图5可见红线所标识的势垒峰的位置上粒子密度不为零。

随着势垒继续向两边移动，粒子继续移向势阱中心，最终又回到如图1中的单个峰的情况，但由于受到两侧势垒的挤压，其峰值要高于图1中的情况，见图6。

随着势垒的继续移动，其对处于系统中心的粒子的影响也逐渐减小（见图7）并最终消失，系统回到最初的纯谐振势作用的情况。

在图8中我们给出随着s 的逐渐增加，高斯势垒的移动，相应系统密度分布变化见图9，从图中可清楚的看到系统的变化。

开始时的势垒处于系统中央区域，粒子分布被劈裂为两个部分。

随着势垒的移动，粒子分布也随之推移。

到达图5中的临界区间时，
粒子所受的谐振
图9. 随着s 的逐渐增加，系统密度分布变化的3D 图。

图10. 系统密度分布变化的临界区域，
两边的粒子分布逐渐回到势阱中央。

势和势垒对能量的贡献达到平衡点，其后粒子回到势阱中央，形成单个峰并最终回到初始的情况。

在图10中给出了临界区域的放大图。

随着势垒的移动，两侧的粒子分布峰逐渐减小，同时在系统中央出现新的峰并逐渐增加，最后两边的粒子分布都回到势阱中央。

五、小结
本文讨论DFT理论的基本原理和方法，从普遍的一位玻色子的DFT出发，讨论了受谐振势和高斯型势垒作用的一维玻色系统两粒子问题，其中高斯势垒由势阱中央出现并成对向系统两侧移动。

对于只有谐振外势作用的情况，系统有精确解，可用以验证计算的有效性[24]。

分布于系统中的粒子受到谐振势和势垒的双重作用时，谐振势和高斯势垒对能量的贡献决定系统的粒子分布。

开始时粒子分布受处于势阱中央的势垒的作用劈裂为两个分布峰，随着高斯势垒的移动，粒子分布随之移动。

当粒子分布移动所引起的谐振势相应能量的增加大于高斯势垒时，粒子开始重新分布于势阱中心并最终全部回到系统中心。

此类处理轻易扩展至更多粒子的情况以及其他有趣的外势作用，同时相互作用的改变也可呈现一些有意思的变化。

对于DFT适用的系统，其计算效率极高，此处的例子正是如此，虽存在高斯型的势垒，但仍得到了不错的结果。

参考文献
[1]L. H.Thomas, Proc. Cambridge Phil. Soc. 23 (5), 542–548 (1927).
[2]Fermi, Enrico, Rend. Accad. Naz. Lincei 6, 602-607 (1927).
[3]P. Hohenberg, W. Kohn, Phys. Rev. 136(3B), B864–B871 (1964).
[4] W. Kohn, L. J. Sham, Phys. Rev. 140, A1133–A1138 (1965).
[5]E. Runge, E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
[6]I. Vasiliev, S. Ogüt, and J. R. Chelikowsky, Phys. Rev. Lett. 86, 1813 (2001).
[7]L. Reining et al., Phys. Rev. Lett. 88, 066404 (2002).
[8]A. Marini et al., Phys. Rev. Lett. 91, 256402 (2003).
[9]B. Walker et al., Phys. Rev. Lett. 96, 113001 (2006).
[10]S. Botti et al., Rep. Prog. Phys. 70, 357 (2007).
[11]D. M. Edwards, J. Magn. Magn. Mater. 45, 151 (1984).
[12]B. L. Györffy and J. B. Staunton, Phys. Rev. Lett. 69, 371 (1992).
[13]A. K. Theophilou, Philos. Mag. B 69, 771 (1994).
[14]V. P. Antropov et al., Phys. Rev. Lett. 75, 729 (1995).
[15]G. M. Stocks et al., Philos. Mag. B, 78, 665 (1998).
[16]S. V. Halilov et al., Phys. Rev. B 58, 293 (1998).
[17]R. Gebauer and S. Baroni, Phys. Rev. B 61, 6459 (2000).
[18]N. D. Mermin, Phys. Rev. 137, A1441 (1965).
[19]N. Argaman , G. Makov, Am. J. Phys. 68, 69 (2000).
[20] N. Argaman, G. Makov, Phys. Rev. B 66, 052413 (2002).
[21] K. Capelle, G. Vignale and B. L. Györffy. Phys. Rev. Lett. 87, 206403 (2001).
[22]E. H. Lieb, W. Liniger, Phys. Rev. 130, 1605 (1963).
[23]Yajiang Hao, Shu Chen, Phys. Rev. A, 80,043608 (2009).
[24]Gao Xianlong et al. Phys. Rev. A,73, 053615 (2006).。