计算材料学：基于DFT的计算晶体弹性模量与等静压计算

您可以点此 下载链接 下载实验报告的 PDF 扫描版本，或者继续阅读本页面，其由 OCR 技术生成但可能存在识别误差的电子文档，建议以扫描版为准，感谢您的查阅。

一、实验预习

1. 实验目的

（1）以MgO晶体为计算对象，学习在Materials Studio（CASTEP模块）中利用弹性常数任务（Elastic Constants Task）直接计算MgO的完整弹性常数矩阵 $C_{ij}$，比较不同超胞类型（Common Cell与Wigner-Seitz Cell）以及不同应变步数（4步、5步、6步）对计算结果的影响，评估计算参数的收敛性与稳定性；

（2）由弹性常数张量出发，通过Voigt-Reuss-Hill（VRH）平均方案推导MgO的多晶弹性参数，包括体弹性模量 $B_0$、杨氏模量 $E$、泊松比 $\nu$、剪切模量 $G$，并与实验值对比分析误差来源；

（3）对弛豫后的MgO平衡结构施加一系列等静压（Isostatic Pressure），在Materials Studio中进行受压几何优化，提取压力$P$与晶胞体积$V$的关系，通过$P$-$V$曲线拟合得到体弹性模量 $B_0$，与弹性常数法和BM方程法的结果进行横向对比，理解三种方法的适用条件与精度差异；

（4）结合等静压概念，理解高压实验（如热等静压HIP工艺）的物理背景，以及DFT计算在指导材料高压合成中的理论意义。

2. 实验原理

2.1 弹性常数的DFT计算——应变能法

对平衡结构施加小应变 $\boldsymbol{\varepsilon}$，体系总能变化为：

$$ \Delta E = E(\boldsymbol{\varepsilon}) - E_0 = \frac{V_0}{2}\sum_{i,j=1}^{6} C_{ij}\varepsilon_i\varepsilon_j + \mathcal{O}(\varepsilon^3) $$

CASTEP弹性常数任务通过自动施加一组线性无关的应变模式，对每种模式取若干步幅（4步、5步或6步），拟合 $\Delta E(\delta)$ 抛物线，提取对应的 $C_{ij}$ 分量。

对于MgO的立方对称性（Fm$\bar{3}$m），独立弹性常数只有三个：

$$ \mathbf{C} = \begin{pmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \end{pmatrix} $$

最大应变幅度设置为0.003，在线性响应范围内保证精度。

2.2 Voigt-Reuss-Hill平均方案

由弹性常数推导多晶弹性参数，分别采用等应变（Voigt上界）和等应力（Reuss下界）假设，取算术平均（Hill平均）：

Voigt上界：

$$ B_V = \frac{C_{11}+2C_{12}}{3}, \qquad G_V = \frac{C_{11}-C_{12}+3C_{44}}{5} $$

Reuss下界（由柔度矩阵 $S_{ij} = C_{ij}^{-1}$ 推导）：

$$ B_R = \frac{1}{3S_{11}+6S_{12}}, \qquad G_R = \frac{5}{4S_{11}-4S_{12}+3S_{44}} $$

其中：

$$ S_{11} = \frac{C_{11}+C_{12}}{(C_{11}-C_{12})(C_{11}+2C_{12})}, \quad S_{12} = \frac{-C_{12}}{(C_{11}-C_{12})(C_{11}+2C_{12})}, \quad S_{44} = \frac{1}{C_{44}} $$

Hill平均：

$$ B_{\text{Hill}} = \frac{B_V+B_R}{2}, \qquad G_{\text{Hill}} = \frac{G_V+G_R}{2} $$

工程弹性参数：

$$ E = \frac{9B_{\text{Hill}}G_{\text{Hill}}}{3B_{\text{Hill}}+G_{\text{Hill}}}, \qquad \nu = \frac{3B_{\text{Hill}}-2G_{\text{Hill}}}{2(3B_{\text{Hill}}+G_{\text{Hill}})} $$

Pugh比（韧脆性判据）：

$$ \frac{B}{G} > 1.75 \Rightarrow \text{韧性材料}; \qquad \frac{B}{G} < 1.75 \Rightarrow \text{脆性材料} $$

2.3 等静压与P-V关系

等静压是指在各方向施加相同大小的静水压力 $P$，体系各向同性压缩。热力学定义：

$$ P = -\frac{\partial E}{\partial V}\bigg|_T = -\frac{1}{3}\text{Tr}(\boldsymbol{\sigma}) $$

在DFT计算中，对弛豫后的平衡结构施加外部压力 $P_{\text{ext}}$，令CASTEP在该压力下进行完全晶胞优化（Cell Optimization: Full），得到受压平衡体积 $V(P)$，构建 $P$-$V$ 数据集。

在小压力范围内，$P$-$V$ 关系满足线性近似：

$$ P \approx B_0 \cdot \frac{\Delta V}{V_0} = B_0 \cdot \frac{V_0 - V}{V_0} $$

即 $P$-$\Delta V/V_0$ 图的斜率即为体弹性模量 $B_0$。更精确的描述采用BM状态方程的压强形式：

$$ P(V) = \frac{3B_0}{2}\left[\left(\frac{V_0}{V}\right)^{7/3} - \left(\frac{V_0}{V}\right)^{5/3}\right]\left\{1 + \frac{3}{4}(B_0'-4)\left[\left(\frac{V_0}{V}\right)^{2/3}-1\right]\right\} $$

此式由BM3 EOS对体积求负导数得到，是高压实验中最常用的状态方程压强形式。

此外，由三次多项式拟合 $E(a)$ 也可间接推导 $B_0$，通过：

$$ B_0 = V_0 \frac{d^2E}{dV^2}\bigg|_{V=V_0} = \frac{1}{9V_0}\left.\frac{d^2E}{da^2}\right|_{a=a_0} \cdot a_0^2 $$

2.4 三种B₀计算方法的比较

二、实验过程

1. 实验过程记录

（1）弹性常数计算——超胞设置

以实验一和实验二中弛豫完成的MgO平衡结构（$a_0 = 4.285$ Å）为初始结构，在CASTEP中选择Elastic Constants任务。分别使用两种超胞构型：

• Common Cell：以原始立方原胞（2原子，$a\times a\times a$）直接施加应变，计算效率高
• Wigner-Seitz Cell（WS Cell）：以Wigner-Seitz超胞（4原子）进行计算，对称性处理更完整

对每种超胞类型，分别测试应变步数 $N = 4, 5, 6$，最大应变幅度统一设为0.003，共计6组计算（截图显示了45Steps和55Steps的CASTEP弹性常数设置界面，应变模式矩阵为 $\{xx:1, yy:1, zz:0, yz:0, xz:1, xy:0\}$ 等）。

计算参数：GGA-PBE泛函，截断能370 eV，K点 $6\times6\times6$，Ultrasoft赝势，Koelling-Harmon相对论修正，不使用体积守恒应变（Use volume-conserving strains 未勾选）。

（2）等静压计算——施加外部压力系列

以弛豫平衡结构为起点，在CASTEP几何优化中通过External pressure/stress设置等静压，压力序列为：

$$ P = 0.5,\ 1.0,\ 1.5,\ 2.0,\ 2.5,\ 3.0,\ 3.5,\ 4.0,\ 4.5,\ 5.0,\ 5.5,\ 6.0,\ 6.5,\ 7.0,\ 7.5,\ 8.0,\ 8.5,\ 9.0,\ 9.5,\ 10.0\ \text{GPa} $$

共17个压力点（加上零压弛豫结构共20个数据点）。每个压力下进行完全晶胞优化（Cell optimization: Full，即允许晶胞形状、体积和原子位置同时弛豫），收敛标准：能量 $< 0.001$ eV/cell，力 $< 100.0$ eV/Å，应力 $< 100.0$ GPa，位移 $< 100.0$ Å，最大迭代步数100步。

从实验截图可见，每个压力点均以独立文件夹（如0.5GPa for all CASTEP GeomOpt）组织，对应.castep输出文件记录收敛后的晶胞体积。由Python脚本自动提取各压力点的 $V(P)$ 数据，构建$P$-$V$数据集并进行三种方式拟合（线性、二次多项式、BM3 EOS）。

三、分析讨论

1. 实验结果与分析

（1）弹性常数计算结果

六组弹性常数计算结果汇总如下（实验数据）：

各弹性参数随步数变化分析：

$C_{44}$ 在所有六组计算中高度一致（$132.99 \sim 133.60$ GPa），标准差不超过0.3 GPa，说明剪切弹性常数对计算参数不敏感，收敛性最好。$C_{11}$ 和 $C_{12}$ 随步数和超胞类型有一定波动（波动约5–15 GPa），是 $B_0$ 分散性的主要来源。Common Cell 系统性高于 WS Cell 约5–10 GPa，可能与两种超胞的应变施加方式和边界条件处理差异有关。

取六组计算的均值作为最佳估计：

$$ \bar{C}_{11} \approx 265\ \text{GPa},\quad \bar{C}_{12} \approx 89\ \text{GPa},\quad \bar{C}_{44} \approx 133\ \text{GPa} $$

$$ \bar{B}_0 = \frac{\bar{C}_{11}+2\bar{C}_{12}}{3} \approx 148\ \text{GPa} $$

韧脆性判断（Pugh比）：

$$ \frac{B}{G}\bigg|_{\text{Common}} \approx \frac{150}{114} = 1.32 < 1.75 \Rightarrow \text{脆性材料} $$

$$ \frac{B}{G}\bigg|_{\text{WS}} \approx \frac{145}{112} = 1.29 < 1.75 \Rightarrow \text{脆性材料} $$

两种超胞均判断MgO为脆性材料，与MgO在室温下为典型离子键陶瓷、解理断裂的实验事实一致。

Born力学稳定性验证：

$$ C_{11} - C_{12} \approx 179\ \text{GPa} > 0 \quad \checkmark $$

$$ C_{11} + 2C_{12} \approx 443\ \text{GPa} > 0 \quad \checkmark $$

$$ C_{44} \approx 133\ \text{GPa} > 0 \quad \checkmark $$

三个Born稳定性条件均满足，MgO岩盐结构在零压下力学稳定。

（2）等静压P-V曲线分析

从$P$-$V$曲线（左）可见，随施加等静压增大，MgO晶胞体积从平衡值 $V_0 \approx 78.61$ Å$^3$ 单调减小，在压力0.5–10.0 GPa范围内体积压缩约3–5%，呈现良好的单调线性趋势，非线性效应较小（压力范围有限）。

三种拟合方法结果：

$$ \text{线性拟合：} B_0 = 162.5\ \text{GPa} \quad (R^2 = 0.9923) $$

$$ \text{二次多项式：} B_0 = 114.3\ \text{GPa} \quad (R^2 = 0.9990) $$

$$ \text{BM3 EOS：} B_0 = 112.2\ \text{GPa} \quad (R^2 = 0.9990) $$

归一化$P$-$\Delta V/V_0$关系（右）斜率即为线性体弹性模量，数据点落在直线两侧且散布较小，说明在本实验的压力范围内（$\leq 10.0$ GPa，体积压缩率 $\leq 6\%$）线性近似基本成立。

从$E(a)$三次多项式拟合中提取的体弹性模量：

$$ B_0^{(\text{poly})} = 177.4\ \text{GPa} $$

三种方法结果汇总对比：

分析：

弹性常数法所得 $B_0 \approx 148$ GPa与实验值160 GPa偏差约7.5%，偏低原因主要是GGA-PBE对MgO弹性常数的系统性低估。等静压线性拟合值162.5 GPa与实验值最接近，但 $R^2 = 0.9823$ 偏低，说明在较高压力下线性假设开始失效，数据出现轻微弯曲。实验一的BM3拟合值（35.9 GPa）严重偏低，原因已在实验一中分析，系扫描范围（3.5–5.0 Å）过宽导致拟合曲率失真，不具可比性。三次多项式估算值177.4 GPa偏高，系多项式拟合对二阶曲率的估计受到远端数据点影响所致。综合来看，弹性常数法是提取 $B_0$ 最为系统且精度可控的方法，等静压法物理图像直观，适合延伸至高压区间的研究。

2. 问题提出与讨论

问题：等静压实验与热等静压（HIP）工艺有何联系？DFT计算对高压材料合成有何指导价值？

实验附图展示了热等静压（Hot Isostatic Pressing，HIP）设备示意图：材料在高压容器中同时承受各向均匀压力（$\leq 2000$ bar $\approx 200$ MPa）和高温（$\leq 2000$ °C），实现致密化烧结或扩散连接。其核心物理过程正是本实验研究的等静压响应——即材料在各向同性压力下的体积压缩与结构响应。

DFT计算的等静压$P$-$V$曲线直接给出：

$$ V(P) \Rightarrow \text{各压力下的平衡晶格常数} \Rightarrow \text{预测X射线衍射峰的压力依赖性} $$

这对高压合成（如金刚石、立方氮化硼、高压氧化物相）具有重要指导意义，具体体现在：

相稳定性预测：通过计算不同压力下各候选相的Gibbs自由能 $G(P,T) = E + PV - TS$，判断相变临界压力：

$$ G_{\alpha}(P^*) = G_{\beta}(P^*) \Rightarrow \text{相变点} $$

体弹性模量与可压缩性：$B_0$ 直接决定材料在工业等静压条件下的体积收缩量，指导HIP工艺参数（压力、温度组合）的选取。

与本实验的对比：本实验压力范围（0.5–10.0 GPa）远超工业HIP（$\leq 0.2$ GPa），更接近超高压合成条件（如金刚石压腔实验GPa至数百GPa量级）。DFT的优势在于可以无损地在任意压力下计算，为实验提供先验指导，尤其对于极端条件下难以进行原位测量的体系。

需要指出的是，本实验DFT计算为零温（$T = 0$ K）结果，未包含声子热压修正。实际HIP工艺在高温下进行，热膨胀和热压效应不可忽略，若要进行定量预测，需在准谐近似（QHA）框架下引入温度修正：

$$ P(V,T) = P_{\text{DFT}}(V) + P_{\text{thermal}}(V,T) $$

其中热压项 $P_{\text{thermal}} = \gamma C_V T / V$，$\gamma$ 为Grüneisen参数，需通过声子计算进一步获得。

作者：GARFIELDTOM
邮箱：coolerxde@gt.ac.cn