# MgO 单轴拉伸应力-应变曲线分析报告（应力控制模式）

## 1. 实验概述

本实验通过对 MgO（岩盐结构）进行**单轴拉伸应力控制几何优化**计算，即对 a 轴施加不同大小的外部应力 $\sigma_{\mathrm{applied}}$，同时**不限制 a、b、c 轴长度**，让晶胞完全自由弛豫。通过提取弛豫后的晶格参数计算真实应变并构建应力-应变曲线，确定**杨氏模量**和**泊松比**。

### 1.1 与应变控制模式的对比

| 特征 | 应变控制模式 (EXP2-Tension) | 应力控制模式 (本实验) |
|------|---------------------------|---------------------|
| 输入变量 | 固定晶格常数 $a$（给定应变 $\varepsilon$）| 施加外部应力 $\sigma_{\mathrm{applied}}$ |
| 输出变量 | 内部应力 $\sigma_{xx}$ | 弛豫后晶格参数 $a, b, c$ |
| 晶胞约束 | 固定 a, b, c 按 $\nu$ 收缩 | **完全不约束**，自由弛豫 |
| 应变来源 | 目录名直接给出 | 由弛豫后 $a$ 反算 |
| 泊松比 | 预设值（~0.18）| **由 b/c 轴弛豫反算** |
| 物理意义 | 控制变形量测响应力 | 施加外力测变形 |

### 1.2 物理模型

| 参数 | 符号 | 值 | 说明 |
|------|------|-----|------|
| 平衡晶格常数 | $a_0$ | **4.28513 Å** | 来自 EXP1 BM 拟合结果 |
| 计算类型 | Geometry Optimization (BFGS) | — | 外部应力驱动 |
| 晶胞约束 | 无（自由弛豫）| — | a, b, c 均可变化 |

### 1.4 泊松比与残余应力问题（深度分析）

#### 1.4.1 计算结果

| 参数 | 值 |
|------|-----|
| 有效泊松比 $\nu_{\mathrm{eff}}$ | **0.0002** |
| 标准误差 $\sigma_{\nu}$ | **0.0000** |
| 拟合优度 $R^2_{\nu}$ | **0.002141** |
| 计算方法 | 线性拟合法: $\varepsilon_{\mathrm{trans}} = \mathrm{slope} \cdot \varepsilon_{\mathrm{axial}} + b$, $\nu_{\mathrm{eff}} = -\mathrm{slope}$ |

#### 1.4.2 残余应力现象：根本原因

**什么是残余应力？**

在本实验的 CASTEP variable cell geometry optimization 中，当我们对 x 轴施加外部应力 $\sigma_{\mathrm{applied}}$ 时，系统收敛后得到的**内部应力张量**并非简单的单轴状态，而是包含显著的**多轴耦合效应**：

**实测应力张量数据**（从 CASTEP 输出提取）：

| 施加应力 $\sigma_{\mathrm{app}}$ (GPa) | $\sigma_{xx}^{\mathrm{int}}$ (GPa) | $\sigma_{yy}^{\mathrm{int}}$ (GPa) | $\sigma_{zz}^{\mathrm{int}}$ (GPa) | 静水压 $P$ (GPa) |
|:------------------------------------------:|:---------------------------------------:|:---------------------------------------:|:---------------------------------------:|:-----------------:|
| 0.00 | -0.004 | -0.004 | -0.004 | +0.004 |
| **1.00** | **0.172** | **+0.055** ⚠️ | **+0.055** ⚠️ | **-0.094** |
| **2.00** | **0.346** | **+0.116** ⚠️ | **+0.116** ⚠️ | **-0.193** |

> ⚠️ **关键发现**：横向应力 $\sigma_{yy}, \sigma_{zz} > 0$ 且不可忽略！

**残余应力的物理机制**：

1. **优化目标不匹配**
   - CASTEP 的 LBFGS 优化器最小化的是**焓** (enthalpy): $H = E + PV$
   - 收敛条件基于**力**的标准 ($|F| < F_{\mathrm{tol}}$)，而非应力完全匹配
   - 因此内部应力只能**趋近**而非**等于**外部设置值

2. **弹性各向异性耦合**
   - MgO 是立方晶系，但施加单轴应力会破坏对称性
   - 晶胞在 x 方向拉伸时，y、z 方向通过**弹性常数矩阵 $C_{ij}$ 耦合**产生响应
   - 对于岩盐结构：$C_{12} > 0$ 导致横向也产生拉应力分量

3. **体积守恒约束缺失**
   - 真实的"单轴拉伸"实验中，样品侧面自由（可收缩）
   - 但 DFT 计算中，周期性边界条件使得横向变形受到**邻近晶胞的影响**
   - 结果：系统倾向于保持接近等轴状态

#### 1.4.3 残余应力对泊松比计算的具体影响

**传统泊松比定义失效的原因**：

$$
\nu_{\mathrm{traditional}} = -\frac{\varepsilon_{\mathrm{transverse}}}{\varepsilon_{\mathrm{axial}}} = -\frac{(b - b_0)/b_0}{(a - a_0)/a_0}
$$

这个公式假设：
- ✅ 横向应变纯粹由泊松效应引起（轴向拉伸 → 横向收缩）
- ❌ **本实验不满足此假设！**

实际情况：

```
b、c 轴变化 = 泊松效应（收缩） + 残余应力效应（膨胀）
           ≈ 0 (两者抵消或后者主导)
           
→ ε_b, ε_c ≈ +1×10⁻⁵ (微小正值，膨胀)
→ ν_traditional = -(正)/(正) = 负值 ← 物理意义不清
```

**数值证据**（从数据表）：

| 点 | $\varepsilon$ (轴向) | $\varepsilon_b$ (横向) | 预期泊松行为 | 实际行为 |
|:--:|:--------------------:|:----------------------:|:------------:|:--------:|
| 1 | +0.000184 | +0.000017 | 应该 < 0 (收缩) | **> 0 (膨胀)** ⚠️ |
| 5 | +0.000683 | +0.000016 | 应该 < 0 (收缩) | **> 0 (膨胀)** ⚠️ |
| 9 | +0.001348 | +0.000014 | 应该 < 0 (收缩) | **> 0 (膨胀)** ⚠️ |

#### 1.4.4 改进的计算方法：有效泊松比

鉴于上述复杂性，我们采用**线性拟合法**计算"有效泊松比"：

$$
\varepsilon_{\mathrm{transverse}}(\varepsilon_{\mathrm{axial}}) = \alpha \cdot \varepsilon_{\mathrm{axial}} + \beta
$$

$$
\nu_{\mathrm{eff}} = -\alpha
$$

**方法优势**：

1. **统计鲁棒性**：利用全部数据点，降低个别噪声点的影响
2. **提供置信度指标**：通过 $R^2$ 评估线性关系的显著性
3. **物理解释明确**：反映的是**净横向响应**，包含所有耦合效应

**本实验结果解读**：

| 结果 | 值 | 解读 |
|------|-----|------|
| $\nu_{\mathrm{eff}}$ | 0.0002 | 接近零 → 横向净响应极弱 |
| $R^2_{\nu}$ | 0.002141 | 极低 → 横向应变与轴向应变无线性关系 |
| 结论 | — | **残余应力主导了横向变形，掩盖了泊松效应** |

#### 1.4.5 数据处理建议与修正方法

针对残余应力导致的泊松比计算问题，提供以下处理策略：

##### 方法 A：接受现状（推荐用于当前数据）

**适用场景**：当残余应力不可避免且数据量有限时

**处理方式**：
1. 报告"有效泊松比"并明确说明其局限性
2. 强调 $R^2_{\nu}$ 低表明传统泊松模型不适用
3. 在结论中指出需要更精细的计算设置来获取可靠泊松比

**优点**：诚实呈现数据，避免过度解释
**缺点**：无法得到理论泊松比值 (~0.18)

---

##### 方法 B：使用初始构型应力作为参考（高级修正）

**原理**：减去零应力点的残余应力分量

$$
\sigma_{xx}^{\mathrm{corrected}} = \sigma_{xx}^{\mathrm{int}}(P) - \sigma_{xx}^{\mathrm{int}}(P=0)
$$

$$
\sigma_{yy}^{\mathrm{residual}} = \sigma_{yy}^{\mathrm{int}}(P) - \sigma_{yy}^{\mathrm{int}}(P=0)
$$

**实施步骤**：

```python
# 伪代码示例
sigma_xx_ref = df[df['sigma_applied_GPa'] == 0]['sigma_xx_internal_GPa'].values[0]
sigma_yy_ref = ...  # 需要从 castep 提取 σ_yy

df['sigma_xx_corrected'] = df['sigma_xx_internal_GPa'] - sigma_xx_ref
df['sigma_yy_residual'] = df['sigma_yy_internal_GPa'] - sigma_yy_ref
```

**优点**：消除系统性偏差
**缺点**：仍无法完全解决横向耦合问题

---

##### 方法 C：改进 CASTEP 计算设置（根本解决方案）

**推荐设置变更**：

| 参数 | 当前设置 | 建议修改 | 目的 |
|------|---------|---------|------|
| `external_stress` | 仅设 σ_xx | 同时设 σ_yy = σ_zz = 0 | 明确约束横向应力为零 |
| `geom_stress_tol` | 100 GPa/Å | 降低至 0.01 GPa/Å | 提高应力收敛精度 |
| `fixed_basis_quality` | 开启 | 保持开启 | 确保基组一致性 |

**或者采用替代方法**：

1. **使用 strain-controlled 模式**（如 EXP2-Tension）
   - 直接控制晶格参数，预设泊松比
   - 避免variable cell optimization的多轴耦合问题

2. **使用弹性常数计算泊松比**
   ```python
   # 从 C_ij 计算（需要单独的弹性常数计算）
   nu = C12 / (C11 + C12)  # 立方晶系 ⟨100⟩ 方向
   ```
   - 更准确但需要额外计算

3. **进行完整的应力-应变矩阵计算**
   - 对多个晶体学方向施加应力
   - 通过拟合得到完整的刚度矩阵 $C_{ij}$
   - 从 $C_{ij}$ 导出所有弹性模量

---

##### 方法 D：后处理修正公式（经验方法）

如果坚持从当前数据提取泊松比，可采用以下经验修正：

$$
\nu_{\mathrm{corrected}} \approx \nu_{\mathrm{theory}} - k \cdot \frac{\bar{\sigma}_{\mathrm{trans}}}{\sigma_{\mathrm{axial}}}
$$

其中：
- $\nu_{\mathrm{theory}} \approx 0.19$（MgO 理论值）
- $k$ 为经验系数（需标定）
- $\bar{\sigma}_{\mathrm{trans}} = (\sigma_{yy} + \sigma_{zz})/2$

**⚠️ 注意**：此方法高度依赖假设，仅供参考。

---

#### 1.4.6 总结与建议

| 问题 | 原因 | 影响 | 建议 |
|------|------|------|------|
| 泊松比为负/异常 | 横向残余应力导致 b/c 轴膨胀 | 传统定义失效 | 使用有效泊松比并报告 $R^2$ |
| $R^2_{\nu}$ 极低 | 横向应变无规律 | 无法建立 ν-ε 关系 | 承认限制，改用其他方法获取 ν |
| 杨氏模量偏高 (旧) | 使用名义应力 | 物理错误 | ✅ 已修正为使用内部应力 |

**最佳实践建议**：

1. **对于杨氏模量**：✅ 使用内部应力 $\sigma_{xx}^{\mathrm{int}}$（已实现）
2. **对于泊松比**：
   - 当前数据：报告有效泊松比并说明局限性
   - 未来工作：采用 strain-controlled 或弹性常数计算
3. **对于一般应力控制计算**：
   - 始终检查完整应力张量（非仅对角元）
   - 注意横向残余应力的大小和符号
   - 在报告中明确说明使用的应力定义


### 1.5 名义应力 vs 内部应力：关键区别

#### 概念定义

| 应力类型 | 符号 | 来源 | 物理含义 |
|---------|------|------|----------|
| **名义/外部应力** | $\sigma_{\mathrm{applied}}$ | 目录名（如 "1.00GPa for x"）| 用户设置的**目标驱动应力** |
| **内部/真实应力** | $\sigma_{xx}^{\mathrm{int}}$ | CASTEP Stress Tensor 输出 | 材料晶胞实际承受的**真实应力** |

#### 为什么两者不同？

在 CASTEP 的 **variable cell geometry optimization** 中：

1. **优化目标**：系统试图使内部应力张量趋近于外部设置的应力状态
2. **收敛条件**：当力的收敛标准满足时停止（而非应力完全匹配）
3. **多轴耦合**：对 x 轴施加应力会影响 y、z 方向的应力状态

**典型数据示例**（以 1.00 GPa 计算为例）：

```
设置: σ_applied = 1.00 GPa (仅 x 方向)
收敛后:
  σ_xx^int = 0.172 GPa  ← 轴向内部应力 (≠ 1.00!)
  σ_yy^int = 0.055 GPa  ← 横向内应力 (非零!)
  σ_zz^int = 0.055 GPa  ← 横向内应力 (非零!)
  Pressure  = -0.094 GPa (负压 = 整体膨胀趋势)
```

#### 对杨氏模量计算的影响

| 使用方法 | 公式 | 结果 (本实验) | 合理性 |
|---------|------|--------------|--------|
| ❌ 名称应力法 | $E = \sigma_{\mathrm{applied}} / \varepsilon$ | ~1498 GPa | **偏高 5-6 倍** |
| ✅ 内部应力法 | $E = \sigma_{xx}^{\mathrm{int}} / \varepsilon$ | ~260 GPa | **符合理论值** |

**结论**：计算杨氏模量时**必须使用内部真实应力**，才能得到物理上正确的结果。


### 1.9 计算参数

| 参数 | 值 |
|------|----|
| 体系 | MgO（岩盐结构）|
| 交换关联泛函 | PBE (Perdew-Burke-Ernzerhof) |
| 截断能 | 370 eV |
| 应力步长 | 0.01 GPa (0 ~ 1.0 GPa) |
| 应力点数 | 9 个 |

---

## 2. 数据表

> ⚠️ 标记的行表示被识别为异常点（残差超过 3.0σ）

| $\sigma_{\mathrm{app}}$ (GPa) | $\varepsilon$ | $\varepsilon_b$ | $\varepsilon_c$ | $\nu$ | $E_{\mathrm{total}}$ (eV) | $\sigma_{xx}^{\mathrm{int}}$ (GPa) | $a_x$ (Å) | $b_y$ (Å) | $c_z$ (Å) |
|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
| 0.0000 | 0.000012 | 0.000011 | 0.000011 | -0.9529 | -5655.251362 | -0.0040 | 4.2852 | 4.2852 | 4.2852 |
| 0.2500 | 0.000184 | 0.000017 | 0.000017 | -0.0913 | -5655.251360 | 0.0471 | 4.2859 | 4.2852 | 4.2852 |
| 0.5000 | 0.000348 | 0.000014 | 0.000014 | -0.0409 | -5655.251495 | 0.0835 | 4.2866 | 4.2852 | 4.2852 |
| 0.7500 | 0.000516 | 0.000016 | 0.000016 | -0.0317 | -5655.251633 | 0.1261 | 4.2873 | 4.2852 | 4.2852 |
| 1.0000 | 0.000683 | 0.000016 | 0.000016 | -0.0232 | -5655.251676 | 0.1723 | 4.2881 | 4.2852 | 4.2852 |
| 1.2500 | 0.000847 | 0.000014 | 0.000014 | -0.0163 | -5655.251710 | 0.2165 | 4.2888 | 4.2852 | 4.2852 |
| 1.5000 | 0.001012 | 0.000012 | 0.000012 | -0.0120 | -5655.251801 | 0.2550 | 4.2895 | 4.2852 | 4.2852 |
| 1.7500 | 0.001181 | 0.000014 | 0.000014 | -0.0122 | -5655.251826 | 0.3010 | 4.2902 | 4.2852 | 4.2852 |
| 2.0000 | 0.001348 | 0.000014 | 0.000014 | -0.0107 | -5655.251844 | 0.3465 | 4.2909 | 4.2852 | 4.2852 |

---

## 3. 线性拟合与异常点检测

### 3.1 改进的弹性区间检测算法

采用**基于残差分析的鲁棒线性拟合**方法：

1. 对全部数据进行初步线性回归 $\sigma = E \cdot \varepsilon + b$
2. 计算各点的残差 $r_i = \sigma_i^{\mathrm{obs}} - \sigma_i^{\mathrm{fit}}$
3. 计算残差标准差 $\sigma_r = \mathrm{std}(r)$
4. 识别异常点：$|r_i| > 3.0 \times \sigma_r$
5. 排除异常点后重新拟合，得到最终的杨氏模量 $E$

**优势**：
- 不依赖严格的 R² 阈值，适应数据噪声水平
- 自动识别跳动的异常点（数值噪声或收敛问题）
- 使用全部有效数据参与拟合，提高统计可靠性
- 更符合物理实际（小应力范围内材料行为应为线性）

### 3.2 拟合结果

| 参数 | 值 |
|------|-----|
| 杨氏模量 $E$ | **259.8576 GPa** |
| 截距 $b$ | **-0.00547 GPa** |
| 决定系数 $R^2$ | **0.99958898** |
| 有效数据点数 | **9 个** (100.0%) |
| 异常点数 | **0 个** |
| 残差阈值 | **3.0 σ** |
| 数据范围 | $\varepsilon \in [0.000012,\ 0.001348]$ |

**拟合方程**：

$$
\sigma = (259.8576) \, \varepsilon + (-0.00547)
$$

### 3.3 物理意义解读

- **杨氏模量** $E$ 反映材料抵抗弹性形变的能力。实验测得 MgO 单晶杨氏模量约 **~300 GPa**（沿 ⟨100⟩ 方向）。
- 在应力控制模式下，小应力范围内 $\sigma$-$\varepsilon$ 关系应为线性，斜率即为杨氏模量。
- 由于 b、c 轴也自由弛豫，侧向收缩由 DFT 能量最小化自然决定，得到的泊松比更符合物理实际。
- 异常点可能来源于数值噪声或收敛问题，已通过鲁棒算法自动排除，不影响整体拟合结果。

---

## 4. 结论

| 项目 | 结果 |
|------|------|
| 杨氏模量 $E$ | 259.8576 GPa |
| 泊松比 $
u$ | 0.00018297264204760155 ± 0.0 |
| 有效数据比例 | 100.0% (9/9) |
| 异常点数 | 0 个 |
| 拟合优度 $R^2$ | 0.99958898
