第一性原理计算方法讲义

第一性原理计算方法 引言 前面讲述的有限元和有限差分等数值计算方法中 求解的过程中需要知道一 些物理参量 如温度场方程中的热传导系数和浓度场方程中的扩散系数等 这 些参量随着材料的不同而改变 需要通过实验或经验来确定 所以这些方法也 叫做经验或者半经验方法 而第一性原理计算方法只需要知道几个基本的物理 参量如电子质量 电子的电量 原子的质量 原子的核电荷数 布朗克常数 波尔半径等 而不需要知道那些经验或半经验的参数 第一性原理计算方法的 理论基础是量子力学 即对体系薛定额方程的求解 量子力学是反映微观粒子运动规律的理论 量子力学的出现 使得人们对 于物质微观结构的认识日益深入 原则上 量子力学完全可以解释原子之间是 如何相互作用从而构成固体的 量子力学在物理 化学 材料 生物以及许多 现代技术中得到了广泛的应用 以量子力学为基础而发展起来的固体物理学 使人们搞清了 为什么物质有半导体 导体 绝缘体的区别 等一系列基本问 题 引发了通讯技术和计算机技术的重大变革 目前 结合高速发展的计算机 技术建立起来的计算材料科学已经在材料设计 物性研究方面发挥着越来越重 要的作用 但是固体是具有 1023数量级粒子的多粒子系统 具体应用量子理论时会 导致物理方程过于复杂以至于无法求解 所以将量子理论应用于固体系统必须 采用一些近似和简化 绝热近似 Born Oppenheimei 近似 将电子的运动和原 子核的运动分开 从而将多粒子系统简化为多电子系统 Hartree Fock 近似将多 电子问题简化为仅与以单电子波函数 分子轨道 为基本变量的单粒子问题 但是其中波函数的行列式表示使得求解需要非常大的计算量 对于研究分子体 系 他可以作为一个很好的出发点 但是不适于研究固态体系 1964 年 Hohenberg 和 Kohn 提出了严格的密度泛函理论 Density Functional Theory DFT 它建立在非均匀电子气理论基础之上 以粒子数密度作为基本变量 r 1965 年 Kohn 和 Sham 提出 Kohn Sham 方程将复杂的多电子问题及其对应的 薛定谔方程转化为相对简单的单电子问题及单电子 Kohn Sham 方程 将精确的 密度泛函理论应用到实际 需要对电子间的交换关联作用进行近似 局域密度 近似 LDA 广义梯度近似 GGA 等的提出 以及以密度泛函理论为基础 的计算方法 赝势方法 全电子线形缀加平面波方法 FLAPW 等 的提出 使得密度泛函理论在化学和固体物理中的电子结构计算取得了广泛的应用 从 而使得固体材料的研究取得长足的进步 第一性原理计算方法的应用 1 体系的能量 进行第一性原理计算前 首先需要确定体系模型 即模型的晶胞和晶胞中 原子的坐标 对于晶体具有周期对称性 具有三个基矢方向和基矢大小 晶格 常数 由于理论计算确定的平衡晶格常数和实验值有所差别 建立模型前需要 确定平衡晶格常数 晶格常数的确定采用如下步骤 通过改变三个基矢的大小 改变单胞的体积 81 119 通过第一性原理 计算可以得到具有不同体积的模型的能量 通过拟合 Murnaghan 方程 得到晶 体的晶格常数以及单胞的能量 0 000 000 00 11 1 B B VVV E VE VB B BVV 其中 为基态平衡体积 为基态下体系的结合能 相对于对应孤 0 V 00 E V 立原子能量 V 为原胞体积 为体模量 为体模量对压强的导数 如课 0 B 0 B 件中图形所示 可以确定在一定体积下体系的能量达到极小值 即体系的基态 能量 所对应的体积为体系的平衡体积 进而可以得到模型三个基矢的大小确 定晶体的平衡晶格常数 这里需要指出的是不同的第一性原理计算方法给出的 能量 代表的物理意义不同 但是本质上都可以反应体系的稳定性 如总能指 构成体系的原子孤立时的能量减去原子成键放出的能量 结合能是以孤立原子 的能量为零点 体系具有的总能 即原子构成晶体时放出的能量 在上面求得的晶格常数的基础上 根据要研究的物理问题 确定体系中包 含原子数目的多少 建立第一性原理计算模型 第一性原理计算的模型通常选 取一个或几个单胞 超单胞 作为模型 选取的模型具有三个基矢方向 应保 证沿着三个基矢方向平移可以构成无限大的晶体 第一性原理计算输入的原子坐标有两种坐标形式 一种是笛卡尔坐标 Cartesian coordinates 一种是分数坐标 fractional coordinates 如对于 Ni3Al 高温合金 具有如图所示的晶体结构 铝原子位于立方体的顶点 镍原 子位于立方体的面心位置 如果取一个单胞作为研究模型 则三个基矢 1 a 分别为 0 0 0 0 0 0 其中为体系的晶格常数 单胞 2 a 3 a aaaa 中包含四个不等价原子 三个 Ni 和一个 Al 采用笛卡尔坐标四个原子的坐标 可表示为 0 0 0 2 2 0 2 0 2 0 2 2 如果采用aaaaaa 分数坐标表示为 0 0 0 1 2 1 2 0 1 2 0 1 2 0 1 2 1 2 迪卡尔 坐标 x y z 和分数坐标 a b c 之间关系为 a b c x y z 其中 a b c 为 1 a 1 a 1 a 一个原子的三个分数坐标 x y z 为该原子的笛卡尔坐标 图中所示各点表示将 晶格常数的大小取不同值时得到的单胞体积作为横轴 而纵轴表示对应体积a 下将原子坐标输入进行第一性原理计算求得的体系能量 拟合后得到 Ni3Al 的 平衡基矢大小以及体系能量 对于研究合金中的掺杂问题 由于掺杂元素的量很少 所以建立的模型需 要取多个单胞 超晶胞 随着模型中原子数目的增加 第一性原理计算方法的 计算量指数增加 对于掺杂量很低的情况 如 0 1 需要模型中至少取 1000 个原子来和实际相符合 这超出了第一性原理计算在目前的计算机上的计算能 力 100 原子左右 所以建立模型时需要考虑能够反映要研究的实际问题就可 以 假设一个超单胞中只存在一个掺杂原子 这样相邻两个超晶胞中掺杂原子 的间距为超晶胞的基矢大小 一般两个原子之间相隔三到四个原子层 原子之 间的相互作用就可以认为非常小了 所以选取八个单胞构成的超单胞就可以基 本反应掺杂量很低的掺杂问题了 将建立好的模型 进行第一性原理计算可以得到体系的总量 对总能进行 变换可以定义体系的内聚能 形成能以及择优占位能 进而可以对掺杂是否有 利于形成 形成掺杂后对体系稳定的影响而进行分析 如课件中所列 内聚能 体系的总能减去所有原子孤立时的能量 即由于原子之间的相互作用而放出能 量 从而内聚能为负值 越小表示形成的体系越稳定 通过和没有掺杂体系的 内聚能相比较可以看出掺杂元素对体系稳定性的影响 形成能 体系的总能减 去体系中各自元素对应的晶体中原子的能量 形成能可表示各种金属组成合金 的能力 另外通过比较掺杂原子替代合金中不同元素原子时体系的形成能可以 得到掺杂原子倾向位于合金的什么位置 这个差值就可以定义为择优占据能 这里需要强调 各种能量是根据要研究的具体问题来定义的 比如我们要就掺 杂原子倾向于位于合金的什么位置 使用总能是不能得到 因为超单胞模型中 各种原子的数目不相同 而每种原子的能量是不一样的 没有可比性 所以定 义了择优占位能 2 电荷密度 电荷密度就是晶体中电子密度的分布 通过电荷密度可以知道晶体中原子 间的成键状态 如金属键 共价键 离子健 van der Waals 键和氢键 为了 更好的表示原子形成晶体后原子间的电荷转移和成键情况 引入差分电荷密度 即两个体系中电荷密度的差值 这两个体系应该具有相同的超单胞 超单胞中 原子类型可以不一样 而原子位置要基本一致 如课件中所示 用 Ni3Al 中一 个 Ni 被掺杂元素替代时体系的电荷密度减去没有掺杂的 Ni3Al 的电荷密度 而 得到差分电荷密度 通过图形可以清楚地看出由于掺杂元素的存在导致的电子 分布状态的改变 再如 我们可以将 Ni3Al 晶体的电荷密度减去 Ni 和 Al 原子 放在超单胞相同位置时孤立 Ni 和 Al 原子的电荷密度 可以得到 Ni3Al 中 Ni 和 Al 原子成键过程中电子密度分布的变化 从而更好的观察原子之间的成键情 况 所以差分电荷密度相当于一个体系的电荷密度相对于另一个体系的电荷密 度改变 目的是为了更好的研究体系中的成键状态 3 能带 能带理论是目前研究固体中电子运动的一个主要理论基础 是在用量子力 学研究金属电导理论的过程中开始发展起来的 最初的成就在于定性地阐明了 晶体中电子运动的普遍性的特点 例如固体为什么会有导体 非导体的区别 晶体中电子的平均自由程为什么会远大于原子的间距等 在半导体技术上 能 带论提供了分析半导体理论问题的基础 有利的推动了半导体技术的发展 能带理论是一个近似的理论 在固体中存在大量的电子 他们的运动是相 互关联的 每个电子的运动都要受其他电子运动的牵连 价电子是人们最关心 的 在原子结合成固体的过程中价电子的运动状态发生了很大的变化 而内层 电子的变化比较小 可以把原子核和内层电子近似看成是一个离子实 价电子 可以看成在一个等效势场中运动 能带理论的出发点是固体中的电子不再束缚 于个别原子 而是在整个固体内运动 称为共有化电子 电子的能量状态从处 于一个电子能级变到在一个能量范围内都会存在 我们下面从对自由电子的能 量讨论 得出能带的表示方法 自由电子的能量 动能 根据德布罗依波长与动量关系 2 2 k EEpm 其中为普朗克常量 为 成为波数或波矢 phkl hhh 2hp2 kp l 所以电子能量可表示为 即电子能量为波矢 k 的函数 如课件中 22 2Ekm h 图形所示 我们讨论固体中电子的能量通常是在 k 空间 倒空间 内进行的 对于晶体具有周期对称性 其对应的倒空间也具有周期对称性 对于一维 的情况 实空间的周期为 a 则倒空间的周期为 定义第一个周期 2 ap 为第一布里渊区 则电子能量在实空间分布的随晶体的周期对 ap ap 称性变化 转化到 k 空间 电子能量随倒空间的周期对称性发生变化 即电子 能量在第一布里渊区随 k 的变化在整个 k 空间中周期性重复 因此可以得到课 件中能量随 k 变化的图形 相当于自由电子的能带经过周期性势场调制后的结 果 所以我们讨论能带只需要考虑电子能量在第一布里渊区随 k 的变化关系就 可以 由于计算所取晶体超单胞形状的不同会导致第一布里渊形状的变化 以 及晶体中原子化合状态的不同导致电子所受到的周期性势场不同 所以电子能 量随 k 的变化关系在不同的晶体中是不同的 即不同的晶体具有不同的能带结 构 从而反映出不同的物理性质 4 电子状态密度 由前面的讨论可知晶体中电子能量状态可以取一定的能量范围 在此能量 范围内在不同的能量区间电子能量状态的多少或填充这些能量状态的电子数目 是不一样的 电子状态密度反映了这个不同 即在一能量区间内电子状态 数 目 的多少 为在 dE 能量区间内电子能量状态的数目 通过对体 ddEhdh 系电子状态密度的分析可以得到晶体中原子间的电子杂化情况 附 第一性原理计算方法的基本原理 第一性原理计算方法的建立是基于对量子力学薛定鄂方程的求解 定态薛 定鄂方程为 其中为粒子的波函 2 2 2 Vx y zEx y z m Y Y h x y zY 数 为波的强度 反映了粒子出现在 x y z 的概率 即概率密度 2 x y zY V 为晶体中粒子受到的势场 E 为体系的能量 为体系的哈密 2 2 H 2 V m h 顿量 即能量算符 例 氢原子的薛定鄂方程及其基态解 考虑到氢原子核外电子受核吸引的势能为 2 0 4 e V rpe 则氢原子体系的哈密顿算符是 22 2 0 24 e H mrpe h 氢原子的薛定鄂方程的直角坐标系的表达式为 22222 222 0 24 e x y zx y zEx y z mxyzr jjj pe h 将薛定鄂方程从直角坐标变换到球坐标 其结果为 22 2 222222 0 1112 sin 0 sinsin4 me rE rxrrrr jjj qj qqqqfpe h 由于氢原子核电场是求对称的 故方程解是球对称的 即 0 j q 0 j f 在这种条件下 氢原子的薛定鄂方程简化为 2 2 22 0 12 0 4 me rE rxrr j j pe h 这是一个系数含变量的二阶齐次微分方程 考虑到电子在无穷远处的几率 为零 以及在整个空间的几率为 1 可以确定方程的解 2 rj 中的 C1 和 C2 得到 为了求解常数 代入简 12 rr C eC e aa j 3 r e a a j p a 化的薛定鄂方程可得 2 2 22 0 222 0 4 mEme rr aa pe hh 所以 2 2 2 0 mE a h 2 2 0 2 20 4 me a pe h 可以得到 为玻尔半径 2 2 00 1 4 me a pea h 0 a 4 22 0 13 6 32 me EeV pe h 0 3 0 1 r e a j pa 通过上面对氢原子中电子的薛定鄂方程的求解 一方面我们了解了如何求 解薛定鄂方程 加深对波函数 本征能量等概念的理解 一方面我们也可以理 解解析求解薛定鄂方程困难 对于氢原子内部只有一个电子 并且上面的例子 中我们只求解了基态能级 对于晶体由电子和原子核组成 电子和原子核都是 运动的 因此体系的能量为电子能量 原子核能量和电子和原子核的相互作用 能 所以哈密顿量可以表示成三部分 如课件中所列 这是每个粒子有三个坐 标 即三个变量 这时对薛定鄂方程进行解析求解显然是不现实的 为了求解 薛定鄂方程 人们提出了一系列的近似和方法 下面将简要论述 1 绝热近似 在固体体系中 由于原子核的质量是电子质量的 103 105倍 所以体系中 电子的运动速度比原子核快得多 可以认为 当核发生一个微小扰动时 迅速 运动的电子可瞬时调整 达到新的平衡 因而在求解电子问题时 可近似认为 原子核固定在给定的位置 这就是所谓的绝热近似 Born Oppenheimer 近似 因此可以认为原子核是固定在给定位置 将原子核的坐标作为参数 而不是变 量 这样可以将薛定鄂方程简化为课件中所列的形式 这样将方程中的变量数 目大大减少 同时体系哈密顿量的形式也得到简化 2 单电子近似 经过绝热近似得到的多电子体系的薛定谔方程中 由于哈密顿量中包含多 体相互作用项 第二项 该项不能分离变量 因而方程难以直接解析求解 为了求解多电子薛定谔方程 需要引入单电子近似 对于含有 N 个电子的多体 系统 假设每个电子都近似的看成是在原子核及其他 N 1 个电子所形成的平均 势场中运动 这样就将多体问题简化成了多个单体问题 若第 i 个电子的波函数为 已包含电子的自旋 则 N 电子体系的波 ii r 函数可由 N 个电子所占据的轨道波函数的乘积来构成 即 Hartree 波函数 化 简后得到单电子满足的薛定鄂方程 方程中为平均势场 即离子实的晶格 i V r 周期势和体系中所有电子产生的平均库仑势 例 分子轨道法计算氢分子的结合能 所谓分子轨道法就是多原子形成分子后 电子不再属于单个原子 也就是说电 子不再在原来的原子轨道上运动 而是在新的分子轨道上运动 以氢分子为例 设想氢分子的氢原子 a 和氢原子 b 当它们是自由原子时 各 自一个价电子 归一化的波函数分别用 a b 氢分子实际是个四体问题 当两个原子相互靠近 波函数交叠 形成氢分子 这时每个电子均为 a 原子和 b 原子共有 哈密顿量为 22 22 12121212 H 22 aabb VVVVV mm hh 其中下脚标 1 2 分别表示两个电子 这里考虑了绝热近似 即原子核认为是静 止的 忽略电子与电子之间的相互左右 上式可以分解为两部分 12 V 22 22 12111222 H H H 22 abab VVVV mm hh 考虑 Hatree 近似 其中 这是 111111 H r E r yy 2 2 111112 1 H 22 ab VVV m h 单电子波动方程 他的解称为分子轨道 11 r y 分子轨道的波函数可以选原子波函数的线性组合 分子轨道波函数应有如下形 式 为归一化系数 r C r r ab yyy r C r r ab yyy CC 由 可以得到体系中电子能量 加上原子核间的相互作用 r H r Edryy 能 即得到分子的总能 减去两个孤立氢原子的能量即得到氢分子的结合能 考虑到电子是 Fermi 子 其波函数应满足反对称条件 即 Pauli 不相容原理 这样就可以将体系的波函数写成 Slater 行列式 化简后得到 Hartree Fock 方程 等式左边最后一项是交换作用势 需要通过自洽求解 Hartree Fock 近似的 重要意义是提出了平均场和单电子近似的概念 在求解过程中利用迭代自洽求 解 这对于以后计算物理的发展起到了深远的影响 但是由于 Hartree Fock 近 似本身仅仅考虑了多体系统中的交换能 而忽略了相关能修正 所以他不能作 为具有相互作用的多电子体系采用单电子近似的严格理论依据 单电子近似的 近代理论基础是在密度泛函理论基础上发展起来 自洽求解流程图如课件中所 示 3 密度泛函理论 1 Hohenberg Kohn 定理 密度泛函理论的基本物理思想是体系的基态物理性质可以仅仅通过电子密 度来确定 由量子力学知道 由哈密顿描述的电子体系的基态能量和基 r H 态波函数都可由能量泛函取最小值来决定 而对于 N EH 电子体系 外部势能完全确定了哈密顿 因此 N 和决定了体系基 V r H V r 态的所有性质 所以 当总粒子数 N 不变时 多电子体系的基态能量是基态密度的唯一泛 函 接下来就是如何对能量泛函作变分处理 并将多体问题严格转化为单E 电子问题 2 Kohn Sham 方程 电子的能量可以表示为 有相互作用粒子系统的动能泛函 外场对电子的作用 电子间库仑排斥作用即 Hartree 项 电子间的交换关联作用 也是的泛函 r Kohn Sham 单电子方程 2 1 2 rKSiii Vrrr 这里 电荷密度用单电子波函数表示 2 ii i rnr 单电子有效势为 KSxc r Vv rdrVr rr Kohn Sham 方程的基本思想就是用无相互作用粒子模型代替有相互作用粒 子哈密顿量中的相应项 而将有相互作用粒子的全部复杂性归入交换关联相互 作用泛函中去 从而导出了单电子 Kohn Sham 方程 Kohn Sham 方程在密度泛 函理论的框架内是严格的 但是交换关联能泛函的具体形式未知 实 xc Er 际中需要对其进行近似 3 局域密度近似 Local Density Approximation LDA 它是处理交换关联泛函的一个简单可行的近似 由 Kohn 和 Sham 提出 目 前得到了非常广泛的应用 其基本想法是 利用均匀电子气的密度 来得 r 到非均匀电子气的交换关联泛函 局域密度近似下 表示为 xc Er LDA xcxc Errdr 其中是密度为的均匀电子气的每个粒子的交换关联能 xc r r LDA 适用于在局域 Fermi 波长和 Thomas Fermi 波长尺度上变化缓 r F TF 慢的体系 80 在实际电子体系中 在局域 Fermi 波长和屏蔽长度的尺度上 电 子密度并不是缓慢变化的 但对于原子 分子和固体的许多基态性质 包括键 长 键角 LDA 计算仍然给出了非常有用的结果 4 超晶胞技术 固体物理中所研究的体系通常具有空间平移对称性 一般可以采用周期性 边界条件的倒易空间能带方法进行处理 第一原理赝势平面波方法 Pseudopotential Plane Wave 就是这类计算方法中的一种 它是非常高效的 ab initio 量子力学计算方法 已经被广泛应用于研究各种材料的晶体结构和电子 结构性质 第一原理赝势方法是以平面波作为基函数 通过构造超原胞 进行 第一原理电子结构计算的方法 1 Bloch 定理 Bloch 定理指出 在具有空间平移对称性的固体体系中 体系的电子波函 数可以写成 ik r ii ref r 这里是与原胞周期性有关的部分 满足如下关系 i f r ii f rf rR 其中 R 是体系的晶格格矢 可以用倒格矢作为波矢的平面波进行展开 i f r iG r ii G G f rce 这里倒格矢由晶体的格矢定义如下 m 为整数 因此 每一G R 2G Rm 个电子波函数都可以用平面波展开 i k G r ii G G rce 在倒格矢取到无穷大时 构成一组完备基 因此单电子波函数按G ik r e 照公式 2 32 可以精确展开 但是在实际工作中 显然不可能将 G 取到