1. 课程导入与分子动力学工程概览
材料在服役中很少只承受单一方向的载荷。桥梁拉索被拉长，混凝土柱被压短，芯片互连铜线在热循环中反复拉伸与压缩，电池极片在充放电过程中经历体积胀缩。宏观上我们看到的是应力—应变曲线、屈服强度、弹性模量；而在原子尺度上，这些现象对应的是键长变化、位错形核、晶界迁移、孔洞长大与相变。分子动力学（Molecular Dynamics, MD）正是把这两端连接起来的计算工具：它用牛顿方程追踪每个原子的运动，再通过统计力学把原子轨迹换算成工程可用的力学量。
本课程的目标不是把 MD 当作黑箱软件来调用，而是把它当作一套“可验证的工程实验”。你需要理解输入文件里每一个参数为什么这样设，理解输出文件里每一条曲线为什么这样走，最终能够独立完成一次拉伸或压缩模拟，并对结果做出可信的工程判断。
一、拉伸压缩模拟的工程意义
为什么工程上需要原子尺度的拉伸压缩模拟？至少有三类问题，连续介质力学给不出满意答案。
第一类是尺寸效应。当构件特征尺寸降到纳米量级，表面原子占比急剧上升，材料的弹性模量、屈服强度、断裂方式都会偏离块体值。纳米线、纳米薄膜、纳米多晶的力学行为必须用原子模型描述。
第二类是缺陷演化机制。宏观塑性理论用位错密度等平均量描述，但位错如何形核、如何与孪晶界和晶界交互、孔洞如何合并，这些机制直接决定材料的韧脆转变。MD 可以逐帧观察这些事件。
第三类是极端条件。高应变率、极低温、强辐照损伤等条件下，实验表征困难，而 MD 可以在受控条件下重复计算，给出趋势和机理。
因此，拉伸压缩 MD 模拟的工程价值可以概括为：为材料设计提供原子尺度的“虚拟力学试验”，为宏观本构模型提供物理参数和机制依据，为实验现象提供机理解释。
二、核心输出量：从原子轨迹到工程曲线
一次拉伸或压缩模拟，最核心的输出是应力—应变曲线。但这条曲线不是直接“算出来”的，而是由若干中间量逐步构造的。理解这条链路，是理解整个课程的关键。
首先是应变。工程上常用工程应变 ε = (L − L₀)/L₀，其中 L₀ 是初始盒长，L 是当前盒长。在 MD 中，我们通常沿某一方向（如 z 轴）以恒定应变率拉伸，每步把盒子长度增加 ΔL = L₀ · ε̇ · Δt，同时让原子按比例仿射变换坐标。这里必须明确参照约定：本课程约定拉伸为正应变、压缩为负应变，应力正负号与之一致，拉应力为正、压应力为负。
其次是应力。MD 中应力通过维里应力公式计算：
σ_αβ = (1/V) [ Σ_i m_i v_iα v_iβ + Σ_i Σ_{j>i} r_ijα f_ijβ ]第一项是动能贡献，第二项是原子间力贡献。对于拉伸压缩模拟，我们最关心的是拉伸方向的正应力分量，例如沿 z 方向拉伸时关注 σ_zz。维里应力对截断半径、系综和温度控制方式敏感，因此后处理时必须说明计算条件。
第三是结构演化量。仅有应力—应变曲线还不够，我们还需要知道材料内部发生了什么。常用的结构分析量包括：
- 径向分布函数 g(r)：判断是否发生非晶化或熔化；
- 共同近邻分析（CNA）与配位数：区分 FCC、BCC、HCP 及无序结构；
- 位错提取算法（DXA）：识别位错线、位错类型与密度；
- 孔洞体积分数与团簇分析：追踪损伤演化；
- 势能—应变曲线：反映键断裂与相变事件。
这些量共同构成“力学响应 + 结构机制”的证据链。只给一条应力—应变曲线，工程判断是不充分的。
三、学习路线：从单晶到多晶，从弹性到断裂
本课程的学习路线遵循“由简到繁、由可控到真实”的原则。建议按以下顺序推进：
- 单晶理想拉伸：建立完美 FCC 单晶，沿 [001] 或 [111] 方向拉伸，观察弹性段、屈服点和位错形核。这是理解 MD 拉伸基本流程的起点。
- 含缺陷单晶：引入空位、间隙原子、预置位错或裂纹，观察缺陷如何影响屈服与断裂。
- 纳米多晶：用 Voronoi 方法构造多晶模型，研究晶粒尺寸对强度的影响（Hall–Petch 与反 Hall–Petch）。
- 压缩与剪切耦合：在压缩条件下观察孪生、相变和孔洞闭合。
- 温度与应变率效应：系统改变温度和应变率，提取激活体积等参数。
- 结果验证与工程换算：把 MD 结果与实验或连续介质模型对比，讨论尺度外推的边界。
每一阶段都要求你能够独立完成建模、加载、计算、后处理和结果判断，而不是只跑通一个示例。
四、工具链：从建模到可视化的完整闭环
MD 拉伸压缩模拟的工具链通常包括四层：建模与前处理、计算引擎、后处理与分析、可视化与验证。
建模层常用 LAMMPS 的 lattice、region、create_atoms 命令构建单晶，用 Atomsk 或 Python 脚本构建多晶和含缺陷结构。计算引擎以 LAMMPS 为主，配合 MEAM、EAM、Tersoff 等势函数描述金属、半导体和共价材料。后处理层用 OVITO 做结构分析和可视化，用 Python（NumPy、Matplotlib、pymatgen）做曲线处理和统计。验证层则把 MD 结果与文献、实验或更高精度方法（DFT）对比。
需要强调的是，工具链的核心不是软件名称，而是数据流：输入结构 → 势函数 → 系综与加载 → 轨迹与热力学输出 → 应力应变与结构分析 → 工程判断。下面这张图给出了这一数据流的整体结构，以及每一层的关键输入输出。
五、从第一性原理到工程判断：本课程的基本要求
在进入具体操作之前，需要建立三条基本认识。
第一，势函数是物理假设，不是万能公式。EAM 适合金属，Tersoff 适合共价键，MEAM 适合方向性键合。选错势函数，后面所有计算都不可信。每次模拟前必须说明势函数来源、截断半径和验证情况。
第二，系综与加载方式决定结果含义。NVT 下拉伸，温度由热浴控制，能量不守恒；NVE 下拉伸，温度会随变形变化。工程上更关心等温条件，因此常用 NVT 或 NPT。应变率通常取 10⁸–10¹⁰ s⁻¹，远高于实验，这是 MD 的固有局限，必须在结果讨论中说明。
第三，后处理与计算同等重要。维里应力的计算对截断半径、原子体积定义和温度修正敏感。同一批轨迹，不同后处理参数可能给出不同的屈服强度。因此，后处理脚本必须可复现、参数必须记录、结果必须交叉验证。