简述第一性计算原理、密度泛函理论和分子动力学

第一性计算原理

第一性原理计算 (First-principles calculations)是从量子力学理论出发，在得知了电子的静止质量 $m_{0}$ ，电子电荷量 $e$ ，普朗克常数 $h$ ，光速 $c$ 以及玻尔兹曼常数 $k_b$ 等基本物理量，便可根据原子核和电子相互作用的原理以及基本运动规律，经过多个近似处理后求解薛定谔方程，进而得到材料几乎所有基态性质的计算方法。不含时的薛定谔方程通常用以下公式表达：

\hat{H}\psi(x)=E\psi(x)

$\hat{H}$ 是哈密顿算符，通常情况下是给定的， $\psi(x)$ 是电子的波函数，从 $\hat{H}$ 出发求解出每个电子所处的轨道波函数 $\psi(x)$ 和相应的轨道能量 $E$ 。求解波函数，需要采用一组已知的函数，这组函数就叫作基组也称为基函数。波函数 $\psi(x)$ 可以写成基函数的线性组合的形式：

\psi(x)=\sum_iX_i(x)c_i

代入薛定谔方程：

\sum_i\hat{H}X_i(x)c_i=E\sum_iX_i(x)c_i

两边同时乘以 $X_j(x)$ 再对 $x$ 积分，把 $\int X_j(x)\hat{H}X_i(x)dx$ 记作 $H_{ji}$ ，采用正交基组：

\sum X_j(x)X_i(x)dx=\delta_{ij}=\begin{cases}1,&i=j\\0,&i\neq j\end{cases} (2-4)

薛定谔方程可化简为：

\sum_iH_{ji}c_i=Ec_j

把 $c_j$ 写成向量的形成，薛定谔方程写成矩阵乘法的形成：

\begin{pmatrix}H_{11}&H_{12}&...\\H_{21}&H_{22}&...\\\vdots&\vdots&\ddots\end{pmatrix}\begin{pmatrix}c_1\\c_2\\\vdots\end{pmatrix}=E\begin{pmatrix}c_1\\c_2\\\vdots\end{pmatrix}

用线性代数可以求解矩阵形式的薛定谔方程，因此对于给定的系统，只要写出它的哈密顿量，再选择一组基组，就可以求解出波函数。

对于含有 n 个电子 N 个原子核的分子，其基态的哈密顿量为：

\hat{H}(r,R)=-\frac{1}{2}\sum_{i}\nabla_{i}^{2}-\frac{1}{2M_{I}}\sum_{I}\nabla_{I}^{2}+\sum_{i<j}\frac{1}{\left|r_{i}-r_{j}\right|}-\sum_{I,i}\frac{Z_{I}}{\left|r_{i}-R_{I}\right|}+\sum_{I<J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}

其中 $M_{\mathrm{l}}$ 为第 $I$ 个原子核的质量， $r_i$ 为电子的坐标， $R_\iota$ 为原子核的坐标， $Z_\iota$ 为第1个原子核的电荷量。 $\nabla_i$ 和 $\nabla_I$ 分别表示对 $r_i$ 和 $R_I$ 求梯度。由于原子核的质量远大于电子的质量，所以电子的速度远大于原子核的速度。研究电子原子核间的运动，可以近似她认为原子核的位置不发生改变，从而将电子的运动和原子核的运动分开处理。这种近似被称为玻恩-奥本海默(Born-Oppenheimer)近似，也称为绝热近似.

虽然绝热近似简化了对薛定谔方程的求解，但对于具有 N 个电子的体系，每个电子有三个空间坐标，用第一性原理的方法描述电子间相互作用过程，需要求解 3N 个维度的薛定谔方程。如果计算机求解薛定谔方程的每个电子坐标取点数为 M,那么对于描述 N 个电子状态的多电子波函数，计算机需要 $M^\mathrm{N}$ 个点数来进行计算和存储数据。当 N>3 时，所需的计算量已远超目前计算机求解薛定谔方程的能力，无法直接求解多电子体系的波函数。

密度泛函理论

1928 年 Hartree 提出除了要考虑的电子外，其余电子间的相互作用平均化这一想法，成功把多电子问题转化为单电子问题。电子通常处于体系的基态即能量最低的量子态，单电子的基态波函数描述为：

\varepsilon_{i}\psi_{i}(r)=\left[-\frac{\hbar}{2m}\nabla^{2}+V_{\mathrm{ext}}(r)+\sum_{j}\int dr^{\prime}\left|\psi_{j}(r^{\prime})\right|^{2}\frac{1}{\left|r_{-}r^{\prime}\right|}\right]\psi_{i}(r)-\sum_{j,i=j}\int dr^{\prime}\frac{\psi_{i}^{*}(r^{\prime})\psi_{i}(r^{\prime})}{\left|r_{-}r^{\prime}\right|}\psi_{j}(r)

其中包含电子与原子核相互作用的势能，电子间库仑势能以及电子间交换作用势能。基于 Hartree-Fock 方程去计算电子结构的经典方法，需要求解复杂的电子波函数。密度泛函理论用电子密度取代了波函数将其作为计算的基本量，把多电子问题转换为对电荷密度及其含时演化的计算[67]。电荷密度仅是三个变量的函数，实际求解处理起来更方便。泛函与函数非常相似，函数的值由一个或者多个变量决定，一个变量只能决定一个函数值。而泛函是函数的函数，其值由一个函数作为变量而决定，且该函数只能决定泛函中的一个取值。1964 年，Hohenberg 和Kohn 等人提出了Hohenberg-Kohn（H-K）第一定理，该定理明确指出了基态波函数和基态电荷密度间存在一一对应的关系，对于同一电子密度，其动能项和相互作用项必然相同[68]。体系的基态能量可以写成以下形式：

E=E[\rho(r)]

其中 $\rho(r)$ 为基态电荷密度。H-K 第二定理指出只要知道泛函的具体形式，就可以通过不断调整电荷密度直到由泛函所确定的能量达到最小化，并找到相应的基态电荷密度。

Kohn 和 Sham 提出了 Kohn-Sham 方程：

\left[-\frac{\nabla^2}{2}+V_{\mathrm{KS}}[\rho(r)]\right]\varphi_i(r)=\varepsilon_i\varphi_i(r)

其中 $V_\mathrm{\kappa s}[\rho(r)]$ 包含了原子核-电子相互作用，Hartree 势(电子与其自身的相互作用)和交换相互作用。Kohn-Sham 方程构造了一个电子间无相互作用的体系，选择合适的泛函可以使得体系的电荷密度分布与真实体系的电荷密度分布相同。基于 Kohn-Sham 方程，多体问题(外部电势产生的电子间相互作用)成功地被化简成一个没有相互作用的电子在有效势场中运动的问题，极大降低了计算量使求解薛定谔方程成为了可能。求解基态问题的 Kohn-Sham 方程是一个自洽的方程组，通过迭代的方法求解。最初给定一个初始电荷密度分布 $\rho_m(r)$ ,一般由原子波函数线性叠加得到。其次根据给定的电荷密度分布依次求出经典库仑势、交换关联势，从而求解Kohn-Sham 方程，继而由 Kohn-Sham 方程求解得出的波函数构建新的电荷密度分布 $\rho_{out}(r)$ ,比较 $\rho_{ln}(r)$ 和 $\rho_{out}(r)$ 是否达到自洽标准，达到则输出结果，未达到则继续下一轮的计算，直到两次电荷密度按分布自治为止。

赝势与交换关联泛函

求解 Kohn-Sham 方程,必须知道对应的 $V_{\mathrm{KS}}[\rho(r)]$ ：

V_{\mathrm{KS}}[\rho(r)]=V_{\mathrm{ext}}(r)+\int\frac{\rho(r')}{\left|r-r'\right|}\mathrm{d}^{3}r'+V_{\mathrm{xc}}[\rho(r)]

对于交换关联项 $V_\mathrm{xc}$ ，交换关联泛函并未给出具体形式，只能用近似的泛函来计算该势能项。常见的有局域密度近似(Local Density Approximation,LDA)和广义梯度近似 (Generalized Gradient Approximation, GGA)。LDA 近似使用均匀电子气来计算体系的交换关联项 $V_\mathrm{xc}$ ，均匀电子气的交换关联能 $E_\mathrm{xc}$ 是可以精确求解的，其余部分采用对自由电子气体进行拟合的方法来处理。GGA 近似是基于对 LDA 近似的修正，引入电荷密度的一阶梯度来解释电荷密度变化较快的部分。 $V_\mathrm{ext}(r)$ 在原子核附近存在数值稳定性的问题，靠近原子核附近，势能会变化得很剧烈，必须采用非常细密的空间网格才能准确描述原子的势能。考虑到原子的最内壳层电子能量被原子核束缚得非常紧，相对外层电子来说可以被认为是不活跃的电子。将内壳层电子看作为屏蔽电子，和原子核一起形成所谓的赝势(Pseudopotential)。赝势的特点是给定一个较小的截断半径 $r_\mathrm{c}$ ，大于该半径的电子波函数及其库仑势是精确的，而小于此半径的电子波函数及其库仑势是平滑的，且不改变电离能。引入赝势后，不需要大量平面波作为基组就能得到较好的结果，大大减少了计算所需的时间和资源。

色散矫正

密度泛函理论在计算固体物理性质等方面得到了广泛的应用，与其他解决量子力学多体问题的方法相比，密度泛函理论给出了令人满意的结果，计算所耗费的时间和资源与实验相比要少许多。尽管密度泛函理论在电子结构计算中相对成熟，但是用来描述分子间相互作用，特别是范德瓦尔斯力，或者计算半导体的能带带隙还是有一定程度上的偏差。对于要考虑范德华作用的体系，可以采用半经验的色散矫正方法（DFT-D）实现，DFT-D 矫正能加到原先 DFT 计算所得的体系能量上就是矫正后的能量。DFT-D 当中最著名的就是 Grimme 提出的 D2 和 D3 校正，DFT-D3 矫正能写为：

E_{disp}^{DFT-D3}=\frac{1}{2}\sum_{A\neq B}\sum_{n=6,8}S_{n}\frac{C_{n}^{AB}}{R_{AB}^{n}}f_{damp,n}(R_{AB})

其中 $R_{AB}$ 表示 A 和 B 两个原子间的距离，n 表示距离的 n 次方， $C$ 是原子间色散矫正系数， $S_n$ 是刻度因子，零阻尼函数 $f$ 的表达形式为：

f_{damp,n}(R_{AB})=\frac{1}{1+6[R_{AB}/(S_{r,n}R_{0}^{AB})]^{-\gamma}}

其中是 $R_0^{AB}$ 原子对的截断半径， $\gamma$ 是预设常数，随着原子间距离减小，阻尼函数逐渐衰减为0，使得 DFT-D3 矫正能在较近距离时精确为 0。

分子动力学

分子动力学（Molecular Dynamics，MD）是从经典物理的统计力学出发，对分子间相互作用势函数及其运动方程进行求解的模拟计算。分析分子运动的行为规律，模拟体系的动力学演化，给出微观不可观测量（分子坐标、速度）与宏观可观测量（体系温度，压强，热容）之间的关系，研究体系的平衡态性质和力学性质，是研究材料内部流体行为、通道运输等现象的有效方法。根据玻恩-奥本海默近似，假设原子核为组成体系的经典粒子，按照经典力学牛顿定律，体系中粒子的受力为：

\vec{F}_{i}(t)=m_{i}\vec{a}_{i}(t)

其中 ${\bar{F}_i}$ 为粒子所受的力， $\bar{a}_i$ 为粒子的加速度，而 $m_i$ 是粒子的质量。力可以表示为势能函数对粒子坐标的一阶导数，加速度可以表示为速度对粒子坐标的一阶导数，写成以下形式：

m_{i}\frac{\partial\vec{\nu}_{i}}{\partial t}=-\frac{\partial U}{\partial\vec{r}_{i}}\\\vec{\nu}_{i}(t)=\frac{\partial\vec{r}_{i}}{\partial t}

其中， $\vec{\nu}_i$ 是粒子的速度， $U$ 为相应的势能函数。这样一个随时间演化的体系，可以根据粒子的坐标去描述其运动轨迹，根据粒子的速度去描述其运动轨迹上的变化，从而完全描述其演化过程。通过给定势函数，赋予体系初始的坐标和速度，可以得到一系列包含整个分子动力学过程的坐标与速度。通过统计粒子的坐标和速度，可以得到体系相应的热力学与动力学性质。

从头算分子动力学（ab initio molecular dynamics, AIMD）为研究电子与原子核相互作用的动力学过程的计算方法。AIMD 是在经验力场分子动力学的基础上发展而来，经验力场分子动力学在其求解方程的过程中，势函数（也称之为“力场”）的选取对输出的结果非常重要。势函数的参数往往是通过拟合实验结果和第一性原理计算结果所得，拟合的性质包括内聚能、晶格常数、弹性模量、点缺陷性质、熔点等等，拟合范围越宽，势函数适用的场景就越广泛。描述简单体系的势函数通常用二体势得出的结果就足够精确，描述稀有气体元素常用Lennard-Jones（L-J）势，金属体系需要考虑自由电子气的贡献，常用嵌入原子势（EAM）等等。但在针对化学环境较为复杂的体系（涉及化学键的断裂），经验力场分子动力学无法准确地描述电子的行为。AIMD 将电子-原子核相互作用体系的动力学问题拆分为电子运动问题和原子核运动问题，对电子的运动用量子力学来描述，对原子核的运动用经典力学来描述。用 DFT 方法对体系进行电子结构的计算，得到基态能量和原子受力，在此基础上拟合势函数去进行分子动力学模拟。常见的系综有 NVT 系综、NVE 系综和 NPT 系综，在本项工作中采用Nose-Hooverthermostat 恒温器去进行分子动力学模拟。