LAMMPS计算原子应力：从Virial公式到应力分布可视化

原子级应力是分子动力学中描述局部力学状态的核心量。与有限元中的连续应力场不同，LAMMPS输出的原子应力是离散的、与体积定义相关的量，正确理解其计算原理和后处理方式是获取可靠应力分布的前提。

Virial应力的定义

LAMMPS中每个原子的应力张量通过Virial公式计算：

σ_ij^α = -(1/V^α) * [m^α * v_i^α * v_j^α + Σ_β r_ij^αβ * f_j^αβ]

其中α标记原子，β遍历所有与α有相互作用的原子。第一项是动能贡献（动量流），第二项是势能贡献（力与距离的叉乘）。V^α是原子α的体积，这是应力计算中最容易出问题的参数——因为原子体积本身没有唯一定义。

LAMMPS默认使用compute stress/atom命令输出原子应力：

compute peratom all stress/atom NULL

输出的是6个分量：σ_xx, σ_yy, σ_zz, σ_xy, σ_xz, σ_yz。需要注意，这个输出值乘以(-1)才是通常意义的应力（LAMMPS内部定义压力为正、应力为负）。另外，输出的值没有除以原子体积，单位是压力×体积（atm*ų），需要手动除以原子体积才能得到真正的应力（atm或bar）。

原子体积的确定方法

由于Virial应力依赖原子体积，而LAMMPS不自动分配体积，用户必须自己选择体积定义方案：

方案一：均分体积

最简单的做法——将总模拟盒子体积除以原子数：

variable vol atom vol/N

适用于均匀体系（如纯金属单相），但在界面、缺陷或非均匀体系中误差很大。

方案二：Voronoi镶嵌

LAMMPS内置compute voronoi/atom（需要VORO++包）：

compute voro all voronoi/atom
variable vol atom c_voro[1]

Voronoi方法为每个原子分配一个多面体体积，物理上更合理，适用于多组分体系和不规则结构。

方案三：径向分布积分

对球形壳层积分径向分布函数g(r)得到原子体积，适用于液态体系：

compute rdf all rdf 100

然后在后处理中积分g(r)到第一极小值，得到配位数和等效原子半径。

应力分布的后处理与可视化

沿z方向的应力剖面

界面或薄膜体系中最常见的分析是应力沿某一方向的分布。步骤如下：

# 按z方向分bin
compute chunks all chunk/atom bin/1d z lower 0.5 units box
fix profile all ave/chunk 100 10 1000 chunks &
    c_peratom[1] c_peratom[2] c_peratom[3] &
    file stress_profile.txt

这会输出每个z-bin内的平均应力分量。注意，输出的仍然是未除体积的值，需要在后处理脚本中除以bin内原子总体积。

von Mises等效应力

对于多轴应力状态，工程上常用von Mises应力判断屈服：

σ_vM = sqrt(0.5 * [(σ_xx-σ_yy)² + (σ_yy-σ_zz)² + (σ_zz-σ_xx)² + 6*(σ_xy²+σ_yz²+σ_xz²)])

这个计算需要在后处理中完成，LAMMPS本身不直接输出von Mises应力。可以用Python脚本从dump文件读取原子应力后计算。

力场类型对原子应力的影响

不同力场的应力计算方式有本质差异：

EAM/MEAM势（金属）：只涉及两体相互作用，Virial计算直接明确。应力精度高，与连续介质力学的对应关系清晰。典型应用：纳米压痕模拟中的应力场分析。

Tersoff/SW势（共价键）：包含三体相互作用，Virial公式需要修改。LAMMPS对多体势的处理方式是将多体力等效为两体力的形式，这在物理上不完全严格，但对于应力分布的定性分析足够。

ReaxFF（反应力场）：力分配复杂，LAMMPS内部将每个键的贡献均分到两个原子。应力计算误差较大，建议用足够大的时间平均来降低噪声。

长程库仑（PPPM）：Ewald求和的长程部分对Virial的贡献需要特殊处理。LAMMPS在kspace_style pppm下会自动加入长程Virial修正，但如果设置了compute stress/atom NULL（不加关键字），长程贡献会被忽略。正确做法是：

compute peratom all stress/atom ke pair bond

不包含NULL，这样LAMMPS会把pair（含长程）的Virial都算进去。

实际案例：纳米压痕中的应力分布

以EAM势描述的Cu纳米压痕模拟为例：

• 样品：20×20×15 nm的Cu块体，底部3层固定
• 压头：球形刚性压头，半径5 nm，以0.1 Å/ps速度下压
• 力场：Cu_EAM (Mishin 2001)
• 时间步长：1 fs（金属单位制）

压头下方1000步的应力统计：

compute peratom all stress/atom NULL
compute voro all voronoi/atom
variable sxx atom -c_peratom[1]/c_voro[1]
variable syy atom -c_peratom[2]/c_voro[1]
variable szz atom -c_peratom[3]/c_voro[1]
dump stress all custom 1000 dump.stress id type x y z v_sxx v_syy v_szz

用OVITO可视化dump文件，按σ_zz着色，可以清楚看到压头正下方出现压缩应力区（正值），周围出现环形拉伸应力区（负值）。最大压缩应力约8-12 GPa，与纳米压痕实验的硬度测量值（6-10 GPa）在同一量级。

统计噪声与时间平均

原子应力的瞬时值涨落极大（典型方差是均值的10-100倍），必须做时间平均才能得到有意义的应力场。经验法则：

• 液态体系：至少10万步（100 ps@1fs步长）的平均
• 固态体系：至少1万步（10 ps）的平均
• 界面/缺陷体系：需要更长时间，因为局部结构弛豫慢

减少噪声的另一种方式是空间平均——将应力投影到粗粒化网格上，每个网格bin内包含50-100个原子，统计误差降低约1/√N倍。