📐 分子动力学入门：铜的杨氏模量计算

2.1   牛顿第二定律 —— MD 的基础方程

MD 模拟的核心就是牛顿第二运动定律：

其中 $\vec{F}_i$ 是第 $i$ 个原子受到的合力，$m_i$ 是其质量，$\vec{r}_i$ 是其位置向量。

而力与势能的关系为：

也就是说，只要知道势能函数 $U$，就可以算出每个原子受到的力，再通过牛顿第二定律算出加速度，最后更新位置和速度。

2.2   势能函数（力场）

势能函数描述了原子之间的相互作用能。常见的有以下几类：

2.2.1   对势（Pair Potential）

最简单的势能形式，只考虑两个原子之间的相互作用：

典型代表是 Lennard-Jones (LJ) 势，常用于稀有气体分子：

其中 $r^{-12}$ 项代表短程排斥（原子核之间的排斥），$r^{-6}$ 项代表长程吸引（范德华色散力）。

2.2.2   嵌入原子法（EAM）

EAM（Embedded Atom Method）是模拟金属最常用的势能函数，其核心思想是：

把每个原子看作嵌入在"电子云海"中的"杂质"。除了原子间的两两作用外，每个原子还会感受到周围所有原子贡献的"背景电子密度"。

其中：

• $F(\rho_i)$ — 嵌入能，原子 $i$ 在电子密度 $\rho_i$ 中的能量
• $\phi(r_{ij})$ — 两体对势
• $f(r_{ij})$ — 电子密度贡献函数

本教程使用的 Cu_Zhou04.eam.alloy 就是 EAM 势的一种实现格式。

2.3   数值积分算法

有了力和加速度，如何更新原子的位置和速度？这就需要数值积分算法。最常用的是 Velocity-Verlet 算法：

具体公式为：

2.4   时间步长的选择

时间步长 $\Delta t$ 是 MD 模拟中非常关键的一个参数：

• 太大会导致数值不稳定，原子"飞出去"
• 太小则模拟太慢，浪费计算资源
• 金属体系通常取 $\Delta t = 1 \text{ fs}$（即 $10^{-15}$ 秒）

经验法则：$\Delta t$ 应小于原子振动周期的 $1/20 \sim 1/30$。

2.5   系综与温度压力控制

MD 模拟中需要控制系统所处的热力学系综。最常用的两种：

2.6   周期性边界条件

MD 模拟只能包含有限数量的原子（通常几千到几百万）。为了让小体系能代表宏观材料，使用周期性边界条件（PBC）：

3.1   什么是杨氏模量？

杨氏模量（Young's Modulus）是材料最基本的弹性常数，描述材料在弹性范围内应力与应变的比值：

其中：

• $\sigma = F/A$ — 工程应力（单位面积上的力）
• $\varepsilon = \Delta L / L_0$ — 工程应变（长度变化率）
• $E$ — 杨氏模量（单位：GPa）

铜的典型杨氏模量：多晶铜 $E \approx 117$ GPa。但单晶铜在不同方向上差异很大——这就是各向异性。

3.2   MD 模拟拉伸法的完整流程

3.3   LAMMPS 输入文件结构

我们将整个模拟分为两个文件：

3.4   关键 LAMMPS 命令解读

构建晶体

# 声明 FCC 晶格，晶格常数 3.615 A
lattice fcc 3.615
# 创建模拟盒子，25x12x12 个晶格重复
region box block 0 25 0 12 0 12 units lattice
create_box 1 box
create_atoms 1 box

设置势函数

# EAM/alloy 格式的势函数
pair_style eam/alloy
pair_coeff * * Cu_Zhou04.eam.alloy Cu

弛豫

# 能量最小化：能量收敛到 1e-10，力的收敛到 1e-12
minimize 1.0e-10 1.0e-12 100000 1000000
# NPT 弛豫：温度 300 K，目标压力 0 bar
fix 2 all npt temp 300.0 300.0 0.1 aniso 0.0 0.0 1.0
run 50000

拉伸

# X 方向恒定应变速率拉伸：1e-4 /ps = 1e8 /s
fix DEFORM all deform 1 x erate 1.0e-4 units box remap x
# Y/Z 方向保持零横向应力
fix NPT_YZ all npt temp 300.0 300.0 0.1 y 0.0 0.0 1.0 z 0.0 0.0 1.0
# 输出应力和应变的值
run 50000

4.1   势函数文件获取

本教程使用 Zhou 2004 的铜 EAM 势函数，这是 NIST 机构维护的高质量势函数之一。

下载步骤

1. 打开 NIST 势函数下载页面
2. 点击页面上的 "Download" 按钮，下载 Cu_Zhou04.eam.alloy
3. 将文件放到工作目录中（和 .lammps 文件同级）

4.2   势函数文件结构解析

EAM/alloy 文件是一个文本文件，由三部分组成：

# Cu_Zhou04.eam.alloy  file structure
# ================================================================
# Line 1: comment (description)
# Line 2: element count, format info, atom types
# Line 3+: data for each element
#   - atomic number, mass, lattice constant, crystal structure
#   - embedding energy function F(rho) table
#   - pair potential phi(r) table
#   - electron density function f(r) table
# After: inter-element pair potential phi(r) table
# ================================================================

4.3   工作目录组织

一个井然有序的工作目录能让你的模拟工作更加高效。推荐结构如下：

Youngs_module/                    (project root)
|
+-- Cu_Zhou04.eam.alloy           (potential file, download needed)
|
+-- 01_build_relax.lammps         (step 1: build + relax)
+-- 02_tensile.lammps             (step 2: tensile test)
+-- 03_analysis.py                (step 3: data analysis)
+-- 04_tensile_110.lammps         (advanced: [110] direction)
+-- run_all.bat                   (one-click run script)
|
+-- doc/                          (documentation)
|   +-- tutorial.html
|
+-- Cu_relaxed.data               (relaxed structure, after run)
+-- stress_strain.dat             (stress-strain data, after run)
+-- tensile.lammpstrj             (trajectory file, after run)
+-- log.relax                     (relax log, after run)
+-- log.tensile                   (tensile log, after run)

4.4   关键参数说明

4.5   LAMMPS 单位制说明

本教程使用 units metal 单位制，这是 LAMMPS 中模拟金属最常用的单位制：

5.1   应力-应变曲线

MD 拉伸模拟完成后，stress_strain.dat 文件中包含每一记录步的应变和应力数据。用 Python 绘制出来，会得到类似下面的曲线：

5.2   线性拟合求杨氏模量

在弹性区域（$\varepsilon \approx 0 \sim 2\%$），应力-应变关系近似为线性：

用最小二乘法拟合这条直线，斜率就是杨氏模量 $E$。

python 03_analysis.py

5.3   LAMMPS 中的应力计算

LAMMPS 使用 Virial 应力公式计算原子级的应力张量：

在 LAMMPS 的 metal 单位制下，pxx 的单位是 bar。转换为 GPa：

负号是因为 LAMMPS 中压力的正方向与工程应力的拉应力方向相反（LAMMPS 中拉应力对应负压）。

5.4   预期结果与误差分析

可能的误差来源

• 应变速率效应：MD 的应变速率 ($10^8$/s) 远高于实验 ($10^{-3}$/s)，可能导致模拟值偏高
• 体系尺寸：体系太小会引入尺寸效应，但 14400 个原子对弹性模量通常已足够
• 势函数精度：不同 EAM 势函数给出的杨氏模量可能相差 5~10%
• 温度波动：300 K 的热振动会使应力-应变曲线有统计噪声

A.1   常见问题 (FAQ)

A.2   推荐阅读

A.3   关键公式速查表