第二章 第一性原理计算方法与软件介绍
19世纪末，科学家们发现经典力学和经典电动力学在描述物质的微观系统时存在明显不足，对实验中的许多现象不能做出真正合理的解释。鉴于此，20世纪初物理学家们在旧量子论的基础上建立了量子力学，主要研究原子、分子、凝聚态物质等内部微观粒子的结构、运动规律等性质，目前已广泛应用于物理、化学、材料等学科领域。随着量子力学理论的不断完善，并结合日趋成熟的计算机技术，量子计算模拟成为了现代科学中必不可少的研究手段之一。第一性原理计算(First-principles calculation)，亦称为从头算(Ab-initio calculation)。该计算方法可根据量子力学基本原理，基于密度泛函理论对材料微观体系的状态和性质进行理论上的预测，且计算过程中不需要使用任何经验参数，只需要一些基本物理量(电子电荷质量e 、电子静止质量m 0、光速c 、普朗克常数h 、波尔兹曼常数k B ) 。本工作所选用的计算程序为Materials Studio 软件中的CASTEP 量子力学模块，该模块是基于密度泛函理论的从头算量子力学程序。本章节将简要的介绍密度泛函理论和CASTEP 计算模块。
2.1密度泛函理论概述
第一性原理主要的研究对象是多原子体系。它依据量子力学原理，且在无任何实验参数引入的情况下，将多原子体系当作由自由电子和原子核组成的多粒子体系进行处理。然而，关于量子力学中多粒子体系处理的出发点则为著名的薛定谔方程(Schrödinger Equation)。Schrödinger 方程是量子力学的一个基本方程，也是第一性原理计算方法的核心，它是由奥地利物理学家薛定谔(Schrödinger) 于1926年提出的。该方程可用于描述微观粒子的运动规律，故亦被称为薛定谔波动方程(Schrödinger Wave Equation)，其定态方程描述如下：
2
[-2μ∇2+V (r )]ψ(r,t ) =i ∂(2-1) ψ(r , t ) ∂t
式中ħ为约化普朗克(Plank)常数；μ和V(r)分别表示粒子质量和势场；r 和t 则为体系中所有电子与原子核的位置坐标；Ψ(r,t)是系统波函数，即运动的微观粒子
在V(r)势场下的波函数。但Schrödinger 方程在描述真实的复杂系统时求解过程非常困难，只能处理氢原子等简单的电子体系。因此，研究者们采用一些方法对Schrödinger 方程进行了科学合理地近似和简化，其中最经典的方法是密度泛函理论(Density Functional Theory，DFT) 。DFT 主要是将多电子的问题转化为单电子进行处理，通过基态电子密度来描述体系的基态物理性质。
2.1.1早期Thomas-Fermi 模型
1927年，托马斯(Thomas)和费米(Fermi)建立了Thomas-Fermi(T-F)模型，该模型提出了一个用电子密度来描述体系动能的表达式。但是，T-F 模型认为电子彼此之间互不影响，忽略了电子-电子的交换相互作用，故其计算精度大大降低，并未得到广泛应用。1930年，Dirac 在T-F 模型上引入了描述电子之间的交换相互作用的局域近似，即Thomas-Fermi-Dirac(T-F-D)理论，该理论表达式描述如下： E TF (n ) =C 1⎰d 3rn (r ) +⎰d 3rV ext (r ) n (r ) +⎰d 3rn (r ) +133n (r ) n (r ' ) d rd r ' 2⎰r -r ' (2-2)
式中第一、二、三、四项分别为系统的电子动能、外场势能、电子-电子交换相互作用能和电子间静电作用能的表达式。T-F-D 理论可以通过最低能量的寻找，来获得基态电子密度。但T-F-D 理论模型仍存在较多的漏洞，并且过于简单化，对于分子、原子和凝聚态物质等微观性质的准确描述还需要进一步完善。
2.1.2 Hohenberg-Kohn定理
对于密度泛函理论，其坚实的理论依据始于Hohenberg-Kohn(H-K)定理[92]。该定理是霍恩伯格(Hohenberg)和科恩(Kohn)于1964年在T-F 模型的基础上提出的，它可以总结归纳为两个定理：定理①证明了体系的基态能量仅仅是电子密度ρ(r)的泛函；定理②证明了以基态密度ρ(r)为变量，将体系能量最小化之后就能够得到基态能量。
根据H-K 定理，系统的基态能量泛函可表示为：
E G [ρ(r )]=E[ρ(r )]+⎰V (r ) ρ(r ) dr (2-3)
其中，第二项代表外场中电子的势能作用；而泛函E[ρ(r )]可表示为：
1ρ(r ) ρ(r ') 'E [ρ(r )]=T[ρ(r )]+⎰⎰+E xc [ρ(r )] (2-4) '2∣r -r ∣
式(2-4)中，第一、二、三项分别代表系统的动能、电子间库仑作用能和电子间交换-关联能。虽然H-K 定理证明了通过基态的电子密度分布函数能够得到体系的总能，但仍不能确定体系总能与基态电子密度分布函数之间所存在的具体泛函形式，并且对于如何利用泛函极值性质求解体系总能的问题尚未得到解决。
2.1.3 Kohn-Sham方程
密度泛函理论真正意义上的广泛应用是通过Kohn-Sham(K-S)方程实现的。由于上述H-K 定理无法确定电子密度函数ρ(r ) 、动能泛函T [ρ(r )]和交换-关联能泛函E xc [ρ(r )]。因此，在1965年科恩(Kohn)和沈呂九(Lu-Jeu Sham)在电子气近似的基础上对H-K 定理做了一些创新性的改进，他们通过证明认为一个多粒子体系的粒子密度数可以由一个简单的单粒子波动方程获得，这对密度泛函理论的实际应用起到了关键作用，该方程即为Kohn-Sham 方程，详细描述如下：
体系的电子密度分布是其组成各电子的波函数的平方和，即：
ρ(r ) =∑|ψi (r ) |2 (2-5)
i =1N
而无相互作用体系的动能泛函T s [ρ(r )]则简单认为是N 个单电子的动能之和：
故可得到： {-ħ2∇+V KS [ρ(r )]}ψi (r )=Eψi i (r ) (2-7) 2m ħ22*T s [ρ(r )]=∑(-∇) ψi (r ) ψi (r ) dr (2-6) 2m i =1N
在DFT 中，体系能量是电子密度ρ(r ) 的泛函：
E [ρ(r )]=T s [ρ(r )]+⎰V ext (r ) ρ(r ) dr +ECoul [ρ(r )]+Exc [ρ(r )] (2-8) 式中E Coul [ρ(r )]为库仑能，而库仑势为： V Coul [ρ(r )]=⎰' (2-9) |r -r |'ρ(r ')
对于H-K 定理中未确定的交换-关联能泛函E xc [ρ(r )]则用电子体系的交换-关联能密度εxc [ρ(r )]来代替。
E xc [ρ(r )]=⎰εxc [ρ(r )]ρ(r ) dr (2-10)
因此，电子交换-关联势为： V xc [ρ(r )]=δExc [ρ(r )]δ[ρ(r )]εxc [ρ(r )] (2-11) =δ[ρ(r )]δ[ρ(r )]
结合式(2-8)，(2-9)和(2-11)可知：
V KS [ρ(r )]=V (r ) +V Coul [ρ(r )]+V xc [ρ(r )] =V (r ) +⎰ρ(r ')
|r -r '|dr '+δExc [ρ(r )] (2-12) δ[ρ(r )]
式(2-5)，(2-7)和(2-12)一起被称为Kohn-Sham 方程。K-S 方程用一个无相互作用的自由电子在有效势场中运动的问题解决了复杂且难处理的多粒子体系问题，该有效势场主要包括外势场和电子间的库仑相互作用，而将有相互作用的粒子的复杂性引入到交换-关联能泛函E xc [ρ(r )]中。因此，寻找准确可靠、易于表达的交换-关联能泛函E xc [ρ(r )]形式是迭代求解K-S 方程的关键。
2.2交换-关联能量泛函
根据Kohn-Sham 方程的描述可知，只有确定准确可靠的交换-关联能泛函E xc [ρ(r )]表达形式，才能使K-S 方程具有实际意义。但交换-关联能E xc [ρ(r )]是电子密度ρ(r ) 的泛函，主要依赖于整个空间的电子密度分布，求解过程较为复杂，迄今为止仍没有准确的表达形式。基于上述问题，研究者们对于E xc [ρ(r )]的处理提出了不同的近似方案，目前较为广泛采用的是局域密度近似(Local Density Approximation ，LDA) 和广义梯度近似(Generalized Gradient Approximation，GGA) 。
2.2.1局域密度近似
局域密度近似(LDA)模型是Kohn 和Sham 在1965年提出的。该近似方法的基本思想是通过均匀电子气密度函数得到非均匀电子气体系的交换-关联能，即将非均匀电子系统的整个空间理想地分为许多足够小的体积元dr ，并认为体积元dr 中电子气体是均匀的无相互作用的。同时，他们考虑了交换-关联能E xc [ρ(r )]与电子自旋程度的关系，故该近似方法也可称为局域自旋密度近似(Local Spin
Density Approximation，LSDA) 。LDA 下的交换-关联能可表示为：
LDA [ρ(r )]=E x LDA [ρ(r )]+E c LDA [ρ(r )]=⎰εx [ρ(r )]ρ(r ) dr +⎰εc [ρ(r )]ρ(r ) dr (2-1 E xc
3)
式中，E x [ρ(r )]和E c [ρ(r )]分别为电子气的交换能和关联能；εx [ρ(r )]和εc [ρ(r )]分别代表交换能密度函数和关联能密度函数；而ρ(r ) 需要考虑电子的自旋，即ρ(r ) 包括自旋向上的电子密度ρ↑(r ) 和自旋向下的电子密度ρ↓(r ) 。
LDA 是一种比较准确的近似方法，它能够对原子、分子和凝聚态物质的许多微观性质给出满意结果，直接推动了密度泛函理论的广泛应用。但LDA 主要适用于均匀电子气或空间中电子密度变化缓慢的体系，对于一些能量梯度较高的情况，该近似方法则存在较大的误差。
2.2.2广义梯度近似
由于LDA 本身存在的局限性，研究者们尝试对其进行了一系列的改进，其中较为成功的修正方案是广义梯度近似(GGA)。GGA 在LDA 的基础上引入了电荷密度的非局域梯度项进行校正，更好地考虑了体系电子密度的不均匀性，对非均匀电子体系的处理具有更高的计算精度。广义梯度近似的表达式描述如下：
GGA E xc [ρ(r )]=⎰εxc [ρ(r ) , |∇ρ(r ) |]ρ(r ) dr (2-14)
其中，∇ρ(r ) 为电荷密度梯度项。
相比于LDA ，GGA 属于半局域化形式，并且在计算精度上有了一定程度的提高，尤其是晶格参数的计算，但其值与实验结果相比仍然存在误差。常用的GGA 形式有PBE(Perdew-Burke-Ernzerhof)、RPBE(Revised-Perdew-Burke-Ernzerhof)、PW91(Perdew-Wang)、WC(Wu-Cohen)等。值得注意的是，采用GGA 进行计算的效果并非总优于LDA ，在具体的科研计算过程中，需要对计算体系进行收敛性测试，通过计算结果与实验数据的对比，选择误差较小、计算效率高的近似方法。
第二章 第一性原理计算方法与软件介绍
19世纪末，科学家们发现经典力学和经典电动力学在描述物质的微观系统时存在明显不足，对实验中的许多现象不能做出真正合理的解释。鉴于此，20世纪初物理学家们在旧量子论的基础上建立了量子力学，主要研究原子、分子、凝聚态物质等内部微观粒子的结构、运动规律等性质，目前已广泛应用于物理、化学、材料等学科领域。随着量子力学理论的不断完善，并结合日趋成熟的计算机技术，量子计算模拟成为了现代科学中必不可少的研究手段之一。第一性原理计算(First-principles calculation)，亦称为从头算(Ab-initio calculation)。该计算方法可根据量子力学基本原理，基于密度泛函理论对材料微观体系的状态和性质进行理论上的预测，且计算过程中不需要使用任何经验参数，只需要一些基本物理量(电子电荷质量e 、电子静止质量m 0、光速c 、普朗克常数h 、波尔兹曼常数k B ) 。本工作所选用的计算程序为Materials Studio 软件中的CASTEP 量子力学模块，该模块是基于密度泛函理论的从头算量子力学程序。本章节将简要的介绍密度泛函理论和CASTEP 计算模块。
2.1密度泛函理论概述
第一性原理主要的研究对象是多原子体系。它依据量子力学原理，且在无任何实验参数引入的情况下，将多原子体系当作由自由电子和原子核组成的多粒子体系进行处理。然而，关于量子力学中多粒子体系处理的出发点则为著名的薛定谔方程(Schrödinger Equation)。Schrödinger 方程是量子力学的一个基本方程，也是第一性原理计算方法的核心，它是由奥地利物理学家薛定谔(Schrödinger) 于1926年提出的。该方程可用于描述微观粒子的运动规律，故亦被称为薛定谔波动方程(Schrödinger Wave Equation)，其定态方程描述如下：
2
[-2μ∇2+V (r )]ψ(r,t ) =i ∂(2-1) ψ(r , t ) ∂t
式中ħ为约化普朗克(Plank)常数；μ和V(r)分别表示粒子质量和势场；r 和t 则为体系中所有电子与原子核的位置坐标；Ψ(r,t)是系统波函数，即运动的微观粒子
在V(r)势场下的波函数。但Schrödinger 方程在描述真实的复杂系统时求解过程非常困难，只能处理氢原子等简单的电子体系。因此，研究者们采用一些方法对Schrödinger 方程进行了科学合理地近似和简化，其中最经典的方法是密度泛函理论(Density Functional Theory，DFT) 。DFT 主要是将多电子的问题转化为单电子进行处理，通过基态电子密度来描述体系的基态物理性质。
2.1.1早期Thomas-Fermi 模型
1927年，托马斯(Thomas)和费米(Fermi)建立了Thomas-Fermi(T-F)模型，该模型提出了一个用电子密度来描述体系动能的表达式。但是，T-F 模型认为电子彼此之间互不影响，忽略了电子-电子的交换相互作用，故其计算精度大大降低，并未得到广泛应用。1930年，Dirac 在T-F 模型上引入了描述电子之间的交换相互作用的局域近似，即Thomas-Fermi-Dirac(T-F-D)理论，该理论表达式描述如下： E TF (n ) =C 1⎰d 3rn (r ) +⎰d 3rV ext (r ) n (r ) +⎰d 3rn (r ) +133n (r ) n (r ' ) d rd r ' 2⎰r -r ' (2-2)
式中第一、二、三、四项分别为系统的电子动能、外场势能、电子-电子交换相互作用能和电子间静电作用能的表达式。T-F-D 理论可以通过最低能量的寻找，来获得基态电子密度。但T-F-D 理论模型仍存在较多的漏洞，并且过于简单化，对于分子、原子和凝聚态物质等微观性质的准确描述还需要进一步完善。
2.1.2 Hohenberg-Kohn定理
对于密度泛函理论，其坚实的理论依据始于Hohenberg-Kohn(H-K)定理[92]。该定理是霍恩伯格(Hohenberg)和科恩(Kohn)于1964年在T-F 模型的基础上提出的，它可以总结归纳为两个定理：定理①证明了体系的基态能量仅仅是电子密度ρ(r)的泛函；定理②证明了以基态密度ρ(r)为变量，将体系能量最小化之后就能够得到基态能量。
根据H-K 定理，系统的基态能量泛函可表示为：
E G [ρ(r )]=E[ρ(r )]+⎰V (r ) ρ(r ) dr (2-3)
其中，第二项代表外场中电子的势能作用；而泛函E[ρ(r )]可表示为：
1ρ(r ) ρ(r ') 'E [ρ(r )]=T[ρ(r )]+⎰⎰+E xc [ρ(r )] (2-4) '2∣r -r ∣
式(2-4)中，第一、二、三项分别代表系统的动能、电子间库仑作用能和电子间交换-关联能。虽然H-K 定理证明了通过基态的电子密度分布函数能够得到体系的总能，但仍不能确定体系总能与基态电子密度分布函数之间所存在的具体泛函形式，并且对于如何利用泛函极值性质求解体系总能的问题尚未得到解决。
2.1.3 Kohn-Sham方程
密度泛函理论真正意义上的广泛应用是通过Kohn-Sham(K-S)方程实现的。由于上述H-K 定理无法确定电子密度函数ρ(r ) 、动能泛函T [ρ(r )]和交换-关联能泛函E xc [ρ(r )]。因此，在1965年科恩(Kohn)和沈呂九(Lu-Jeu Sham)在电子气近似的基础上对H-K 定理做了一些创新性的改进，他们通过证明认为一个多粒子体系的粒子密度数可以由一个简单的单粒子波动方程获得，这对密度泛函理论的实际应用起到了关键作用，该方程即为Kohn-Sham 方程，详细描述如下：
体系的电子密度分布是其组成各电子的波函数的平方和，即：
ρ(r ) =∑|ψi (r ) |2 (2-5)
i =1N
而无相互作用体系的动能泛函T s [ρ(r )]则简单认为是N 个单电子的动能之和：
故可得到： {-ħ2∇+V KS [ρ(r )]}ψi (r )=Eψi i (r ) (2-7) 2m ħ22*T s [ρ(r )]=∑(-∇) ψi (r ) ψi (r ) dr (2-6) 2m i =1N
在DFT 中，体系能量是电子密度ρ(r ) 的泛函：
E [ρ(r )]=T s [ρ(r )]+⎰V ext (r ) ρ(r ) dr +ECoul [ρ(r )]+Exc [ρ(r )] (2-8) 式中E Coul [ρ(r )]为库仑能，而库仑势为： V Coul [ρ(r )]=⎰' (2-9) |r -r |'ρ(r ')
对于H-K 定理中未确定的交换-关联能泛函E xc [ρ(r )]则用电子体系的交换-关联能密度εxc [ρ(r )]来代替。
E xc [ρ(r )]=⎰εxc [ρ(r )]ρ(r ) dr (2-10)
因此，电子交换-关联势为： V xc [ρ(r )]=δExc [ρ(r )]δ[ρ(r )]εxc [ρ(r )] (2-11) =δ[ρ(r )]δ[ρ(r )]
结合式(2-8)，(2-9)和(2-11)可知：
V KS [ρ(r )]=V (r ) +V Coul [ρ(r )]+V xc [ρ(r )] =V (r ) +⎰ρ(r ')
|r -r '|dr '+δExc [ρ(r )] (2-12) δ[ρ(r )]
式(2-5)，(2-7)和(2-12)一起被称为Kohn-Sham 方程。K-S 方程用一个无相互作用的自由电子在有效势场中运动的问题解决了复杂且难处理的多粒子体系问题，该有效势场主要包括外势场和电子间的库仑相互作用，而将有相互作用的粒子的复杂性引入到交换-关联能泛函E xc [ρ(r )]中。因此，寻找准确可靠、易于表达的交换-关联能泛函E xc [ρ(r )]形式是迭代求解K-S 方程的关键。
2.2交换-关联能量泛函
根据Kohn-Sham 方程的描述可知，只有确定准确可靠的交换-关联能泛函E xc [ρ(r )]表达形式，才能使K-S 方程具有实际意义。但交换-关联能E xc [ρ(r )]是电子密度ρ(r ) 的泛函，主要依赖于整个空间的电子密度分布，求解过程较为复杂，迄今为止仍没有准确的表达形式。基于上述问题，研究者们对于E xc [ρ(r )]的处理提出了不同的近似方案，目前较为广泛采用的是局域密度近似(Local Density Approximation ，LDA) 和广义梯度近似(Generalized Gradient Approximation，GGA) 。
2.2.1局域密度近似
局域密度近似(LDA)模型是Kohn 和Sham 在1965年提出的。该近似方法的基本思想是通过均匀电子气密度函数得到非均匀电子气体系的交换-关联能，即将非均匀电子系统的整个空间理想地分为许多足够小的体积元dr ，并认为体积元dr 中电子气体是均匀的无相互作用的。同时，他们考虑了交换-关联能E xc [ρ(r )]与电子自旋程度的关系，故该近似方法也可称为局域自旋密度近似(Local Spin
Density Approximation，LSDA) 。LDA 下的交换-关联能可表示为：
LDA [ρ(r )]=E x LDA [ρ(r )]+E c LDA [ρ(r )]=⎰εx [ρ(r )]ρ(r ) dr +⎰εc [ρ(r )]ρ(r ) dr (2-1 E xc
3)
式中，E x [ρ(r )]和E c [ρ(r )]分别为电子气的交换能和关联能；εx [ρ(r )]和εc [ρ(r )]分别代表交换能密度函数和关联能密度函数；而ρ(r ) 需要考虑电子的自旋，即ρ(r ) 包括自旋向上的电子密度ρ↑(r ) 和自旋向下的电子密度ρ↓(r ) 。
LDA 是一种比较准确的近似方法，它能够对原子、分子和凝聚态物质的许多微观性质给出满意结果，直接推动了密度泛函理论的广泛应用。但LDA 主要适用于均匀电子气或空间中电子密度变化缓慢的体系，对于一些能量梯度较高的情况，该近似方法则存在较大的误差。
2.2.2广义梯度近似
由于LDA 本身存在的局限性，研究者们尝试对其进行了一系列的改进，其中较为成功的修正方案是广义梯度近似(GGA)。GGA 在LDA 的基础上引入了电荷密度的非局域梯度项进行校正，更好地考虑了体系电子密度的不均匀性，对非均匀电子体系的处理具有更高的计算精度。广义梯度近似的表达式描述如下：
GGA E xc [ρ(r )]=⎰εxc [ρ(r ) , |∇ρ(r ) |]ρ(r ) dr (2-14)
其中，∇ρ(r ) 为电荷密度梯度项。
相比于LDA ，GGA 属于半局域化形式，并且在计算精度上有了一定程度的提高，尤其是晶格参数的计算，但其值与实验结果相比仍然存在误差。常用的GGA 形式有PBE(Perdew-Burke-Ernzerhof)、RPBE(Revised-Perdew-Burke-Ernzerhof)、PW91(Perdew-Wang)、WC(Wu-Cohen)等。值得注意的是，采用GGA 进行计算的效果并非总优于LDA ，在具体的科研计算过程中，需要对计算体系进行收敛性测试，通过计算结果与实验数据的对比，选择误差较小、计算效率高的近似方法。