您可以点此 下载链接 下载实验报告的 PDF 扫描版本，或者继续阅读本页面，其由 OCR 技术生成但可能存在识别误差的电子文档，建议以扫描版为准，感谢您的查阅。
一、实验预习
1. 实验目的
（1）掌握基于密度泛函理论（DFT）的第一性原理计算方法，以MgO晶体为计算对象，学习在Materials Studio（CASTEP模块）中进行截断能收敛测试、K点收敛测试的完整流程；
（2）通过对不同晶格常数下MgO体系总能的系统计算，获得 $E$-$a$ 和 $E$-$V$ 数据集，分别采用三次多项式拟合与Birch-Murnaghan（BM）三阶状态方程拟合，提取MgO的平衡晶格常数 $a_0$ 和体弹性模量 $B_0$；
（3）理解BM状态方程的物理意义及其与三次多项式拟合的本质区别，培养计算结果的分析与误差讨论能力。
2. 实验原理
2.1 密度泛函理论基础
DFT将多体电子问题映射为等效的单电子Kohn-Sham方程：
$$ \left[-\frac{\hbar^2}{2m}\nabla^2 + V_{\text{ext}}(\mathbf{r}) + V_{\text{Hartree}}(\mathbf{r}) + V_{\text{xc}}(\mathbf{r})\right]\psi_i(\mathbf{r}) = \varepsilon_i\psi_i(\mathbf{r}) $$
其中交换关联势 $V_{\text{xc}}$ 采用广义梯度近似（GGA-PBE）处理。电子波函数用平面波基组展开，截断条件为：
$$ \frac{\hbar^2}{2m}|\mathbf{k}+\mathbf{G}|^2 \leq E_{\text{cut}} $$
布里渊区积分采用Monkhorst-Pack方案离散采样：
$$ E_{\text{total}} \approx \sum_{\mathbf{k}_i} w_i f(\mathbf{k}_i) $$
2.2 收敛性测试
计算前必须对截断能 $E_{\text{cut}}$ 和K点网格分别做收敛测试，判据为相邻点每原子能量差：
$$ |\Delta E^{\text{atom}}| = \left|\frac{E_i - E_{i-1}}{N_{\text{atom}}}\right| \times 10^3 < \varepsilon, \quad \varepsilon = 1\ \text{meV/atom} $$
两者须按顺序进行：先固定较大K点网格扫描截断能，再固定收敛截断能扫描K点，逻辑依据是截断能决定波函数质量，波函数质量影响每个K点的能量准确性。
2.3 能量最低原理与晶格常数确定
根据量子力学变分原理，体系在平衡结构下总能最低。预先设定一组晶格常数 $\{a_i\}$，对每个 $a_i$ 做固定体积单点能计算，得到一组 $(a_i, E_i)$ 或 $(V_i, E_i)$ 数据，通过拟合找到能量极小值对应的 $a_0$。
2.4 三次多项式拟合
对 $E(a)$ 进行三次多项式拟合：
$$ E(a) = Aa^3 + Ba^2 + Ca + D $$
极小值条件 $\frac{dE}{da}=0$：
$$ 3Aa_0^2 + 2Ba_0 + C = 0 $$
联立可得 $Ca_0^4 + 2Ba_0^2 + 3A = 0$，求解得平衡晶格常数 $a_0$，对应平衡能量 $E_0 = E(a_0)$。
2.5 Birch-Murnaghan三阶状态方程
BM3 EOS是描述固体在压缩/膨胀下能量与体积关系的标准状态方程：
$$ E(V) = E_0 + \frac{9V_0B_0}{16}\left\{(\eta-1)^3 B_0' + (\eta-1)^2(6-4\eta)\right\} $$
其中：
$$ \eta = \left(\frac{V_0}{V}\right)^{2/3} $$
拟合参数及其物理含义：
| 参数 | 含义 |
|---|---|
| $E_0$ | 平衡态总能（eV/f.u.） |
| $V_0$ | 平衡原胞体积（Å$^3$/f.u.） |
| $B_0$ | 体弹性模量（GPa） |
| $B_0'$ | $B$ 对压强的一阶导数（无量纲） |
体弹性模量定义为：
$$ B_0 = -V\frac{\partial P}{\partial V}\bigg|_{V=V_0} = V_0\frac{\partial^2 E}{\partial V^2}\bigg|_{V=V_0} $$
对于立方结构MgO（岩盐型，Fm$\bar{3}$m），$V_0 = a_0^3/4$（每个化学式单元）。
二、实验过程
1. 实验过程记录
（1）建立MgO晶体结构模型
在Materials Studio中建立MgO岩盐结构（空间群Fm$\bar{3}$m），初始晶格常数参考实验值 $a = 4.211$ Å，原胞含2个原子（1个Mg，1个O），Wyckoff位置分别为4a $(0,0,0)$ 和4b $(0.5,0.5,0.5)$。计算设置：GGA-PBE泛函，Ultrasoft赝势，Koelling-Harmon相对论修正，Gaussian展宽0.1 eV。
（2）截断能收敛测试
固定K点网格为 $3\times3\times3$，在310–375 eV范围内每隔5 eV做一次单点能计算，共34个点。提取每原子总能，计算相邻点差值 $|\Delta E^{\text{atom}}|$，以1 meV/atom为收敛阈值。
从截断能收敛曲线可见，总能随截断能增大单调下降后趋于平稳；从差值图（图2）可见，355→360 eV区间首次满足 $|\Delta E^{\text{atom}}| < 1$ meV/atom，此后保持收敛。综合考虑精度与计算代价，确定截断能为360 eV。
（3）K点收敛测试
固定截断能370 eV（略大于收敛值以保证余量），依次测试 $3\times3\times3$、$4\times4\times4$、$5\times5\times5$、$6\times6\times6$、$7\times7\times7$、$8\times8\times8$、$9\times9\times9$、$10\times10\times10$ 共8组K点网格。
从K点收敛曲线可见，$3\times3\times3$ 时能量偏高，$4\times4\times4$ 起能量迅速下降并趋于稳定，$6\times6\times6$ 之后能量变化极小。考虑到MgO为宽带隙绝缘体，费米面简单，K点收敛较快，确定K点网格为 $6\times6\times6$（不可约K点28个）。
（4）BM方程晶格常数扫描
以收敛参数（截断能370 eV，K点 $6\times6\times6$）对一组晶格常数 $a = 3.5, 3.6, 3.7, \ldots, 5.0$ Å（共16个点）分别做单点能计算，获得 $(a_i, E_i)$ 数据集，同时换算为 $(V_i, E_i)$（其中 $V = a^3/4$ per f.u.）。
三、分析讨论
1. 实验结果与分析
截断能收敛测试结果：
MgO截断能扫描结果显示，总能在310 eV时约为 $-1413.53$ eV，随截断能增大单调下降，在360 eV附近进入平台区，最终趋向约 $-1413.69$ eV。差值图中，首次满足 $|\Delta E^{\text{atom}}| < 1$ meV/atom的区间为355→360 eV，首次收敛截断能确认为360 eV，中位数 $|\Delta E^{\text{atom}}| = 1.754$ meV/atom，扫描总体呈单调收敛特征，无异常跳点，SCF收敛质量良好。
K点收敛测试结果：
K点从 $3\times3\times3$（4个不可约K点）扫描至 $10\times10\times10$（110个不可约K点），能量在 $4\times4\times4$ 时发生大幅跳变（从约 $-1413.69$ eV降至约 $-1413.80$ eV），此后趋于稳定。这一跳变反映了 $3\times3\times3$ 采样过稀、布里渊区积分严重不足的情况。$6\times6\times6$ 至 $10\times10\times10$ 间能量变化在0.01 eV量级以内，满足收敛要求。MgO作为绝缘体，K点收敛速度快，$6\times6\times6$ 即可满足精度需求。
三次多项式拟合结果：
对16个晶格常数点进行三次多项式拟合：
$$ E(a) = -3.80691\,a^3 + 54.26085\,a^2 + (-255.30572)\,a + (-1016.65580) $$
$$ R^2 = 0.99887917 $$
由极值条件 $3Aa_0^2 + 2Ba_0 + C = 0$ 求解得：
$$ \boxed{a_0 = 4.284\ \text{Å},\quad E_0 = -1413.865\ \text{eV}} $$
BM3 EOS拟合结果：
将数据转换为每化学式单元的 $(V, E)$ 形式，代入BM3状态方程拟合：
$$ E_0 = -353.4528\ \text{eV/f.u.}, \quad V_0 = 19.876\ \text{Å}^3/\text{f.u.} $$
$$ \boxed{B_0 = 35.89\ \text{GPa},\quad B_0' = 3.9639} $$
$$ R^2 = 0.99995305 $$
由 $V_0$ 换算晶格常数：$a_0 = (4V_0)^{1/3} = (4\times19.876)^{1/3} \approx 4.284$ Å，与三次多项式拟合结果一致。
与实验值对比：
| 参数 | 本次计算（GGA-PBE） | 实验值 | 偏差 |
|---|---|---|---|
| $a_0$ (Å) | 4.284 | 4.211 | +1.73% |
| $B_0$ (GPa) | 35.89 | 160–165 | 偏低 |
晶格常数偏大约1.73%，符合GGA系统性高估晶格常数的规律（GGA低估结合能 $\Rightarrow$ 键长偏长）。体弹性模量偏低较多，主要原因在于：本次扫描的晶格常数范围（3.5–5.0 Å）过宽，包含了大量远离平衡点的非物理区间，导致BM拟合所用数据点的能量曲率信息被稀释；若将扫描范围收窄至平衡值 $\pm5\%$ 以内（即约4.0–4.5 Å），BM拟合精度将显著提升，$B_0$ 将更接近实验值。此外，GGA泛函本身对体弹性模量存在系统性低估。
2. 问题提出与讨论
问题：为何三次多项式拟合与BM3 EOS拟合得到的晶格常数一致，但体弹性模量差异较大？
三次多项式拟合的极值位置 $a_0$ 由一阶导数为零决定，对数据的整体趋势敏感，因此两种方法的 $a_0$ 结果一致性较好。但体弹性模量 $B_0$ 正比于 $\frac{\partial^2 E}{\partial V^2}|_{V_0}$，对极小值附近曲线的局部曲率高度敏感。三次多项式是纯经验拟合，对曲率的描述精度不如BM3 EOS这类有物理约束的状态方程。当扫描范围过宽时，极小值两侧的远端数据点会拉动三次多项式的曲率估计，造成 $B_0$ 失真。严格提取 $B_0$ 应采用BM3 EOS + 收窄至平衡值$\pm5\%$范围内的扫描数据，这是与实验值可比的标准做法。
关于真空层退化问题的延伸讨论（见实验截图中的二维材料示意图）：报告附图展示了两层二维材料（Layer 1蓝色，Layer 2红色）在真空层厚度 $c = 10.66$ Å下的波函数密度（Overlap factor = 0.0049）。该overlap factor虽然已经很小，但仍不为零，说明10.66 Å的真空层尚未完全消除层间波函数耦合。计算表明，Interaction strength约为 $\sim 0$ meV，已进入近似退耦（Decoupled）区间，但对于高精度功函数或激发态计算仍需增大至15–20 Å以确保物理隔离。
作者：GARFIELDTOM
邮箱：coolerxde@gt.ac.cn