LAMMPS 安装与入门：从选安装路线到跑通第一个拉伸模拟

第 1 步

安装路线怎么选

先按“要不要 MPI 并行”“要不要 CUDA”“有没有 root”三个条件选路线。选错的典型后果是装完才发现不能并行，或者 GPU 根本没用上。

版本

用 30Sep2026 还是 22Jul2025：先看编译器

2026 年 9 月 30 日发布的新稳定版提高了编译工具链门槛。集群上的 GCC 或 CMake 偏旧时，直接编新版会在 configure 阶段失败。

判断方法：先运行 gcc --version 和 cmake --version。GCC 低于 10.4 又要用 KOKKOS，就选 stable_22Jul2025_update6，或者先装新编译器。CMake 版本不够时，可以用 python3 -m pip install --user cmake 装一个新版，不需要 root。

如果 in 文件要和导师或已发表工作的结果对比，记录 lmp -h 输出的完整版本号。30Sep2026 修正了 ARM64 和开启 FMA 的平台上与平台相关的结果差异，两个版本的数值可能有末位差别。

Windows

Windows：安装包够用到什么程度，什么时候换 WSL2

LAMMPS-GUI 安装包适合学习：编辑器能直接运行当前 in 文件并显示 thermo 曲线和快照。它的局限有三条：不能 MPI 并行；GPU 只有 OpenCL 版 GPU 包，需要显卡驱动带 OpenCL 运行时；KOKKOS 只能用串行或 OpenMP。安装包是自签名的，浏览器和 SmartScreen 会提示风险，在提示中选择“仍要运行”即可。

Windows 下多核运行有两种写法。用 GUI 安装包时用 OpenMP：lmp -sf omp -pk omp 8 -in in.file。用 MSMPI 安装包时用 MPI：mpiexec -n 8 lmp -in in.file。官方说明两者叠加通常更慢，二选一即可。

要 CUDA、要编自己的包、或者要和集群保持同一套脚本，就换 WSL2。CUDA on WSL 只在 Windows 端装 NVIDIA 驱动（R495 及以后），WSL 内只装 cuda-toolkit-12-x 元包。不要装 cuda、cuda-12-x 或 cuda-drivers 元包，它们会覆盖 WSL 映射进来的驱动。WSL2 的 GPU 支持要求 Pascal 及以后的显卡。

# Windows PowerShell（管理员）
wsl --install -d Ubuntu-24.04
wsl --set-default-version 2

# 进入 Ubuntu 后：只装 CUDA Toolkit，不装 Linux 显卡驱动
# 按 NVIDIA 下载页选择 Linux > x86_64 > WSL-Ubuntu，安装 cuda-toolkit-12-x 元包
nvidia-smi                 # 若找不到：/usr/lib/wsl/lib/nvidia-smi

# 工作目录放在 Linux 文件系统里，不要放 /mnt/c
mkdir -p ~/md && cd ~/md

第 2 步

Linux / WSL2 源码 CMake 编译（CPU 版，非 root）

basic.cmake 预设打开 KSPACE、MANYBODY、MOLECULE、RIGID 四个包，30Sep2026 还会打开 GRAPHICS。EAM 在 MANYBODY 里，本页的拉伸算例只需要这个包。

# 依赖（Ubuntu 22.04/24.04 或 WSL2 Ubuntu）
sudo apt update
sudo apt install -y build-essential cmake git libopenmpi-dev openmpi-bin python3-pip
cmake --version            # 30Sep2026 要求 >= 3.20；不够就：python3 -m pip install --user cmake

# 获取源码：stable 分支 = 最新稳定版（2026-10 时为 30Sep2026）
git clone -b stable --depth 1 https://github.com/lammps/lammps.git lammps
# 需要与 2025.07.22 一致时改用：git clone -b stable_22Jul2025_update6 --depth 1 ...

cd lammps && mkdir build && cd build
cmake -C ../cmake/presets/basic.cmake \
      -D BUILD_MPI=yes -D BUILD_OMP=yes -D PKG_OPENMP=yes \
      -D CMAKE_INSTALL_PREFIX=$HOME/.local \
      ../cmake
cmake --build . -j 8
cmake --install .            # 装到 ~/.local，不需要 root

# 写进 ~/.bashrc
export PATH=$HOME/.local/bin:$PATH
export LAMMPS_POTENTIALS=$HOME/.local/share/lammps/potentials

# 验证：必须看到 MPI 库名、MANYBODY 等包
lmp -h | sed -n '/MPI v/p;/Installed packages/,/List of individual/p'

• lmp -h 的输出里要有 MPI 库名。CMake 找不到 MPI 开发包（libopenmpi-dev）时会静默编出串行版。这时 mpirun -np 4 实际跑 4 个互不相干的串行副本：日志显示 1 by 1 by 1 MPI processor grid，每行输出重复 4 次。
• mpirun 必须来自编译时用的那套 MPI。集群上 module load 的 MPI 和编译时的不一致时，同样会出现上面的 N 个副本现象。
• 额外的包用 -D PKG_<NAME>=yes 打开，例如 -D PKG_REAXFF=yes、-D PKG_MEAM=yes、-D PKG_EXTRA-COMPUTE=yes。需要大部分常用包时用 -C ../cmake/presets/most.cmake。
• 改了包的组合后，在同一个 build 目录重新执行 cmake 再编译即可，已有设置会保存在 CMakeCache.txt 里。

第 3 步

GPU 编译：GPU 包和 KOKKOS 选哪个

两者都能在 NVIDIA 卡上跑 EAM，但分工方式不同：GPU 包把 pair 计算放到 GPU，其余步骤留在 CPU；KOKKOS 尽量让整个时间步留在 GPU 上。

GPU_API 的默认值是 opencl。只写 -D PKG_GPU=on 会编出 OpenCL 版，在 NVIDIA 卡上也能跑，但不是 CUDA 后端。这一项是网上教程最常漏写的。

GPU_ARCH 和 Kokkos_ARCH 按显卡计算能力填写：V100 是 sm_70 / VOLTA70；T4 和 RTX 20 系列是 sm_75 / TURING75；A100 是 sm_80 / AMPERE80；RTX 30 系列和 A40 是 sm_86 / AMPERE86；RTX 40 系列和 L40S 是 sm_89 / ADA89（不是 sm_90）；H100 是 sm_90 / HOPPER90；RTX 50 系列是 sm_120 / BLACKWELL120。CUDA 13 已不能为 Maxwell、Pascal、Volta 离线编译。22Jul2025 的默认值是 sm_50，配 CUDA 13 时必须显式指定 GPU_ARCH。

KOKKOS 编出的程序换到主版本号不同的 GPU 架构上（例如 7.x 换到 8.x）会直接报错退出，同主版本只会在启动时多一次 JIT 编译。为多种卡准备程序时，按最旧的那张卡编译，或者为每种卡各编一个。

# 先查显卡计算能力：8.6 -> sm_86，8.9 -> sm_89
nvidia-smi --query-gpu=name,compute_cap --format=csv

# 方案 A：GPU 包（CUDA 后端）
mkdir build-gpu && cd build-gpu
cmake -C ../cmake/presets/basic.cmake \
      -D PKG_GPU=on -D GPU_API=cuda -D GPU_ARCH=sm_86 -D GPU_PREC=mixed \
      -D BUILD_MPI=yes -D BUILD_OMP=yes -D PKG_OPENMP=yes \
      -D CMAKE_INSTALL_PREFIX=$HOME/.local \
      ../cmake
cmake --build . -j 16

# 方案 B：KOKKOS（CUDA 后端）；两个 -C 预设可以叠加
mkdir build-kokkos && cd build-kokkos
cmake -C ../cmake/presets/basic.cmake \
      -C ../cmake/presets/kokkos-cuda.cmake \
      -D Kokkos_ARCH_AMPERE86=yes \
      -D BUILD_MPI=yes \
      ../cmake
cmake --build . -j 16

• 经验做法（MatSci 论坛多帖）：用 CMake 时不要再在 src 目录执行 make yes-gpu 或单独编译 lib/gpu。两种构建方式混用时 CMake 会报错，处理办法是在 src 目录执行 make no-all purge 或重新解压源码，只用 -D PKG_GPU=on。
• 经验做法：运行时报 CUDA driver version is insufficient 一类错误，多数是 CUDA Toolkit 比驱动新。nvidia-smi 右上角显示的是驱动支持的最高 CUDA 版本，nvcc --version 显示的是编译用的 Toolkit 版本，前者应不低于后者。
• 编译时找不到 nvcc，把 CUDA 的 bin 目录加入 PATH（例如 export PATH=/usr/local/cuda/bin:$PATH）后，删掉 build 目录重新 cmake。

第 4 步

运行命令：多核、GPU 包与 KOKKOS 的写法

# 纯 CPU，8 个 MPI 进程
mpirun -np 8 lmp -in in.cu_tensile

# 纯 CPU，OPENMP 包，1 个进程 8 线程（Windows GUI 安装包也可用）
lmp -sf omp -pk omp 8 -in in.cu_tensile

# GPU 包：4 个 MPI 进程共用 1 张卡
mpirun -np 4 lmp -sf gpu -pk gpu 1 -in in.cu_tensile

# KOKKOS：1 张卡配 1 个 MPI 进程
lmp -k on g 1 -sf kk -in in.cu_tensile
# Pascal 及以后的卡可以试半邻居表，常更快
lmp -k on g 1 -sf kk -pk kokkos newton on neigh half -in in.cu_tensile
# 多卡且 MPI 不支持 GPU-aware 时（发行版自带 OpenMPI 常见）
mpirun -np 2 lmp -k on g 2 -sf kk -pk kokkos gpu/aware off -in in.cu_tensile

• 用 -in 文件名，不用 < 重定向。部分 MPI 实现下 < 在并行时读不到输入。
• conda 的 mpich 变体在环境的 bin 目录里自带 mpirun（MPICH Hydra），激活环境后直接用 mpirun -np 8 lmp -in in.cu_tensile，不要混用系统里的 OpenMPI。lmp -h 的 MPI 一行会显示 MPICH Version。本页在 macOS 上实测：开着 VPN 时 MPICH 会选中 utun 虚拟网卡，报 OFI poll failed (default nic=utun...)，先 export FI_PROVIDER=tcp 即可正常运行。
• conda 的 CPU 包不含 GPU 和 KOKKOS 包。本页实测：加 -k on 报 Cannot use -kokkos on without KOKKOS installed，加 -sf gpu -pk gpu 1 报 Package gpu command without GPU package installed；-sf omp -pk omp N 可以使用。
• -sf gpu 会给所有有 /gpu 版本的命令加后缀，eam/alloy 会变成 eam/alloy/gpu；-sf kk 同理。不想改 in 文件时，用命令行开关即可切换 CPU 与 GPU。
• KOKKOS 默认假设 MPI 支持 GPU-aware。单卡单进程不受影响。多进程出现段错误，或提示 Turning off GPU-aware MPI since it is not detected 时，加 -pk kokkos gpu/aware off。
• GPU 包的混合精度（mixed）是默认值。发表前用 double 精度编译或 CPU 版跑一段同样的模拟，对比能量和应力。
• 估算耗时：先把 run 改成 2000 步，看日志末尾 Performance 一行给出的 timesteps/s，再按总步数外推。

第 5 步

第一个 in 文件：单晶铜 [100] 单轴拉伸

体系是 32000 个原子的周期性块体铜，用 Mishin 2001 EAM 势。流程是：能量最小化 → 300 K 零压平衡 20 ps → x 方向以 1×10⁹ s⁻¹ 拉伸 300 ps，同时 y、z 方向保持零压。

本页实测（LAMMPS 22Jul2025 update6，conda-forge CPU 版，macOS arm64）：把 region 改成 0 10 0 10 0 10（4000 个原子），其余参数不变，完整跑完 20 ps 平衡和 300 ps 拉伸。300 K 平衡后 lx 为 36.31 Å，对应晶格常数 3.631 Å。plot_ss.py 给出 0–1% 应变的模量 67.3 GPa，与 [100] 杨氏模量的理论值一致：用同一势在 0 K、0.1% 应变下算得 C11 = 169.7 GPa、C12 = 122.3 GPa（Mishin 2001 论文给出 169.9 和 122.6 GPa），由 E = (C11−C12)(C11+2C12)/(C11+C12) 得 67.2 GPa。第一次屈服应力 8.80 GPa，对应应变 0.113，随后应力突降到约 4.3 GPa；σyy、σzz 始终在 ±0.3 GPa 内，说明 y、z 方向的零压控制有效。864 个原子（6×6×6 晶胞）得到 68.6 GPa 和 8.85 GPa（应变 0.115）。这个屈服应力是无缺陷周期性晶体的均匀形核应力，比宏观实验值高约两个数量级，属于理想晶体模型的正常结果。

单核耗时（同一版本与机器，4000 个原子，各跑 1000 步，CPU 占用 99.9%）：NPT 平衡段 204 步/秒（0.82 M atom-step/s），拉伸段带 fix deform、fix print 和 dump，为 156 步/秒，即每步 6.4 ms（0.62 M atom-step/s）。按此外推，32000 个原子单核每步约 51 ms，平衡加拉伸共 32 万步约 4.5 小时。验证期间本机同时运行其他任务，多进程并行的加速比受干扰，本页不给出。

# in.cu_tensile —— 单晶铜 [100] 单轴拉伸
# 适用：LAMMPS 22Jul2025 与 30Sep2026；需要 MANYBODY 包
# ---------- 1. 初始化 ----------
units           metal          # 长度 Å，时间 ps，能量 eV，压力 bar
dimension       3
boundary        p p p
atom_style      atomic

variable        T     equal 300.0               # 温度 K
variable        erate equal 1.0e9*1.0e-12       # 1e9 s^-1 换算成 1/ps

# ---------- 2. 建模 ----------
lattice         fcc 3.615                       # 与势文件给出的 Cu 晶格常数一致
region          box block 0 20 0 20 0 20        # 20x20x20 个晶胞，32000 个原子
create_box      1 box
create_atoms    1 box

# ---------- 3. 势函数 ----------
pair_style      eam/alloy
pair_coeff      * * Cu_mishin1.eam.alloy Cu     # 文件在 LAMMPS 的 potentials/ 目录
neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

# ---------- 4. 能量最小化，同时把盒子压力弛豫到 0 ----------
thermo          100
thermo_style    custom step pe press lx ly lz
fix             relax all box/relax iso 0.0 vmax 0.001
min_style       cg
minimize        1.0e-12 1.0e-12 10000 100000
unfix           relax

# ---------- 5. 300 K、0 bar 平衡 20 ps ----------
reset_timestep  0
timestep        0.001                           # 1 fs
velocity        all create ${T} 4928459 mom yes rot yes dist gaussian
fix             eq all npt temp ${T} ${T} 0.1 iso 0.0 0.0 1.0
thermo          1000
thermo_style    custom step temp pe press lx ly lz
run             20000
unfix           eq

# ---------- 6. x 方向恒应变速率拉伸，y、z 方向保持 0 bar ----------
variable        tmp equal lx
variable        L0 equal ${tmp}                 # 立即求值，记下拉伸前的长度
reset_timestep  0
fix             nptyz all npt temp ${T} ${T} 0.1 y 0.0 0.0 1.0 z 0.0 0.0 1.0 drag 1.0
fix             pull all deform 1 x erate ${erate} units box remap x

variable        strain equal (lx-v_L0)/v_L0
variable        sxx equal -pxx/10000            # bar 转 GPa；取负号使拉应力为正
variable        syy equal -pyy/10000
variable        szz equal -pzz/10000
fix             out all print 100 "${strain} ${sxx} ${syy} ${szz}" file stress_strain.txt screen no title "# strain sxx_GPa syy_GPa szz_GPa"

compute         peatom all pe/atom
dump            traj all custom 2000 tensile.lammpstrj id type x y z c_peatom
thermo_style    custom step v_strain temp v_sxx v_syy v_szz pe press
restart         50000 tensile.*.restart
run             300000                          # 300 ps，工程应变达到 0.30
write_data      final.data

1. 01

   units metal

   决定后面所有数字的单位：距离 Å，时间 ps，能量 eV，压力 bar。EAM、Tersoff、SW 等势文件大多按 metal 单位给出。

2. 02

   lattice fcc 3.615 与 region/create_box/create_atoms

   lattice 同时设定晶格和长度单位，region 里的 0 20 指 20 个晶胞，约 72 Å。晶格常数用势函数的平衡值，否则一开始就有残余应力。Cu_mishin1.eam.alloy 文件头写的是 3.615 Å。

3. 03

   pair_style eam/alloy 与 pair_coeff * * 文件 Cu

   eam/alloy 读 setfl 格式（.eam.alloy），只能写一条 pair_coeff * *，后面按原子类型顺序写元素名，元素名的个数等于原子类型数。setfl 文件含质量，所以不需要 mass 命令。.eam 文件（funcfl 格式）要配 pair_style eam。两者混用会读错文件。

4. 04

   neighbor 与 neigh_modify

   metal 单位默认 skin 为 2.0 Å。自 2Aug2023 起 neigh_modify 默认 delay 0 every 1 check yes。旧教程里的 delay 10 可能产生 Dangerous builds 警告，建议删掉。

5. 05

   fix box/relax + minimize

   最小化时同时调整盒子尺寸，使 0 K 压力为 0。vmax 0.001 限制每步体积变化，防止最小化发散。

6. 06

   velocity + fix npt iso

   在 300 K、0 bar 下平衡 20 ps，让晶格热膨胀到位。跳过这一步直接拉伸，曲线起点会带热应力。阻尼参数 0.1 ps（Tdamp）和 1.0 ps（Pdamp）分别是 100 步和 1000 步，符合手册建议的量级。

7. 07

   variable L0 equal ${tmp}

   ${} 会在读到这一行时立即求值，所以 L0 固定为拉伸前的长度。直接写 variable L0 equal lx 的话，L0 会随盒子一起变，应变恒为 0。

8. 08

   fix npt y/z + fix deform x erate

   deform 以恒定工程应变速率拉长 x 方向。erate 的单位是 1/时间单位，metal 下是 1/ps，所以 1×10⁹ s⁻¹ 要写 0.001。y、z 用 npt 控制为 0 bar，允许泊松收缩，得到的是单轴应力状态。只写 deform 而 y、z 固定，得到的是单轴应变，应力偏高。拉伸方向 x 不要同时用 npt 控压，否则控压会抵消拉伸；这是 MatSci 论坛上拉伸曲线异常最常见的原因之一。

9. 09

   variable sxx equal -pxx/10000

   pxx 是压力张量分量，单位为 bar，拉伸时为负。取负号使拉应力为正，除以 10000 换成 GPa。变量里用到 pxx、pyy、pzz 时，thermo_style 里必须有 press 一类的压力关键字，LAMMPS 才会计算压力张量。本页实测：去掉 thermo_style 末尾的 press，拉伸段的 run 一开始就报 ERROR: Thermo keyword pxx in variable requires thermo to use/init press。

10. 10

    fix print、dump、restart

    fix print 的字符串用双引号包住，里面的 ${strain} 会在每次输出时求值。dump 输出到 .lammpstrj，OVITO 和 VMD 都能直接打开。restart 每 50000 步写一次续跑文件。

变体

改成纳米线：几何、约束与应力换算

纳米线拉伸是块体脚本最常见的改法。需要改三处：建模、控压方向，以及应力归一化体积。

y、z 方向保留周期边界，但留出远大于截断半径（Mishin 势为 5.5 Å）的真空层，纳米线就不会和自己的周期镜像相互作用。上例盒宽约 72 Å，线径约 40 Å，真空约 32 Å。

pxx 是全盒子的维里应力除以盒子体积。盒子体积包含真空，直接用 -pxx/10000 会明显低估纳米线的应力。上面的换算假设截面在拉伸中不变。也可以对 compute stress/atom 求和，再除以原子数乘以单原子体积。无论用哪种，都要在论文中写明归一化体积的定义。

本页实测（同一版本；为缩短时间把长度改为 20 个晶胞，共 7540 个原子；300 K 平衡 4 ps 后以 1×10⁹ s⁻¹ 拉伸到应变 0.03）：平衡后 lx 为 71.52 Å，比建模长度 72.30 Å 短 1.1%，这是表面应力使纳米线轴向收缩的结果，所以 L0 要在平衡之后取。0–1% 应变拟合的模量，用上面的体积换算为 55.8 GPa，直接用 -pxx/10000 只有 13.2 GPa，两者之比约 4.2，等于盒子截面积与纳米线截面积之比。对 compute stress/atom 求和，除以原子数乘单原子体积 a³/4（并按轴向伸长修正），得到的应力与上面的换算相差小于 0.5%。stress/atom 的符号与压力相反，求和后不要再取负号。

# 把第 2 步建模替换为下面几行：[100] 铜纳米线，直径约 4 nm，长约 20 nm
lattice         fcc 3.615
region          box block 0 55 -10 10 -10 10    # 晶格单位；x 取整数个晶胞
create_box      1 box
region          wire cylinder x 0 0 5.5 INF INF # 半径 5.5 个晶格常数，约 19.9 Å
create_atoms    1 region wire

# 第 4 步改为只在 x 方向弛豫盒子
fix             relax all box/relax x 0.0 vmax 0.001
# 第 5 步改为只在 x 方向控压
fix             eq all npt temp ${T} ${T} 0.1 x 0.0 0.0 1.0
# 第 6 步：y、z 是真空，不能再控压，改用 NVT
fix             nvt1 all nvt temp ${T} ${T} 0.1
fix             pull all deform 1 x erate ${erate} units box remap x

# pxx 按整个盒子体积（含真空）归一，要换算到纳米线自身体积
variable        R    equal 5.5*3.615
variable        sxx  equal -pxx*vol/(PI*v_R^2*lx)/10000   # GPa，忽略泊松收缩

势函数

势函数怎么选、从哪里下载

Cu 的 EAM 势不需要另外下载：Cu_mishin1.eam.alloy（Mishin 等，PRB 63, 224106, 2001）、Cu_u3.eam（Foiles 1986）和 Cu_zhou.eam.alloy 都在 LAMMPS 的 potentials/ 目录里。NIST 原子间势库（ctcms.nist.gov/potentials/system/Cu）上，同一个 Mishin 势的文件名是 Cu01.eam.alloy。NIST 对 Mendelev 2013 的 Cu2 势注明了改进堆垛层错能，对 Zhou 2004 的 ipr2 版本注明了修正表格函数的伪振荡。研究位错和孪晶时，按这些说明选择势函数，并在方法部分写明来源。

用 conda 安装时没有这个目录，从 GitHub 上与版本对应的 tag 下载即可：curl -LO https://raw.githubusercontent.com/lammps/lammps/stable_22Jul2025_update6/potentials/Cu_mishin1.eam.alloy，放到 in 文件所在目录，或放到一个固定目录后把 LAMMPS_POTENTIALS 指向它。

OpenKIM 的势不用下载文件，在 in 文件里直接引用模型 ID：kim init EAM_Dynamo_MishinMehlPapaconstantopoulos_2001_Cu__MO_346334655118_006 metal，接着在建模之后写 kim interactions Cu。前提是编译了 KIM 包。conda 版和 apt 版自带 KIM，Windows 预编译包不带。

新版会读取势文件第一行的 UNITS 标签。in 文件的 units 与之不符时，报错 Potential file ... requires metal units but real units are in use。处理方法是修改 units，或换用对应单位的势文件。

单位

metal 与 real 单位，以及最常见的换算错误

• 从 metal 换到 real 时，erate 数值要再乘 10⁻³：1×10⁹ s⁻¹ 在 real 下是 1.0e-6（1/fs）。
• compute stress/atom 的单位是 压力×体积（bar·Å³），并且符号与压力相反。它不是单原子应力，除以单原子体积后才能当应力看。
• thermo 中的 press 是 (pxx+pyy+pzz)/3，单轴拉伸时它约为 -σxx/3。不能用 press 代替 σxx。
• MD 拉伸常用的应变速率在 10⁷–10¹⁰ s⁻¹ 之间，比实验高出很多个数量级，屈服应力会偏高。经验做法：正式计算前在 10⁹ 和 10⁸ s⁻¹ 各跑一次，看屈服应力差多少，再决定是否值得用更低的速率。比较不同尺寸或温度时，保持应变速率一致，并在文中说明。

第 6 步

后处理：thermo 输出、OVITO 与应力-应变曲线

log.lammps 中每个 run 的 thermo 段以 Step 开头的表头行开始，以 Loop time of 一行结束。表头列名就是 thermo_style 里的关键字，变量会显示为 v_strain、v_sxx。Loop time 之后有时间分解表（Pair、Neigh、Comm、Output 等的占比）：Pair 占比低而 Comm 或 Output 占比高时，说明进程数过多或输出过于频繁。

OVITO 能直接打开 tensile.lammpstrj（也能读 .gz）。在 Pipeline 中添加 Common neighbor analysis 修改器，能按 FCC、HCP、Other 给原子着色。在 FCC 铜中，单层 HCP 原子对应孪晶界，两层相邻的 HCP 原子对应本征堆垛层错。VMD 用 vmd -lammpstrj tensile.lammpstrj 打开，它读不了压缩文件。

本页实测：用上面的 4000 原子算例运行 plot_ss.py，输出 E(0-1%) = 67.3 GPa 和 first yield = 8.80 GPa at strain 0.113。小晶胞在第一次屈服后会重新加载，864 原子算例在应变 0.23 处出现 9.9 GPa 的更高峰，所以脚本找的是第一次明显回落之前的最大值，不能直接取全局最大值。lammps.formats.LogFile 读出的 runs 依次对应最小化、平衡和拉伸三段，runs[-1] 的键为 Step、v_strain、Temp、v_sxx、v_syy、v_szz、PotEng、Press。

位错密度随应变的曲线是中文社区常见的后处理需求，CSDN 上有多篇文章演示，脚本多在 QQ 群里分发。下面的 dxa_density.py 只用 OVITO Python 模块（pip install ovito）的公开接口：DislocationAnalysisModifier 输出位错线总长 DislocationAnalysis.total_line_length（Å）和盒子体积 DislocationAnalysis.cell_volume（ų），两者相除再乘 10²⁰ 得到 m⁻²；DislocationAnalysis.length.1/6<112> 是 Shockley 不全位错的长度。

本页实测（OVITO 3.16.1，上面 4000 原子算例的 151 帧拉伸轨迹，用时约 11 s）：应变 0.116 之前 DXA 没有找到位错，FCC 原子比例在应变 0.100 时为 0.936，0.116 时降到 0.625；应变 0.118 时第一次检出位错，密度 4.1×10¹⁷ m⁻²，HCP 原子比例升到 0.149；此后密度在 0 到 1.4×10¹⁸ m⁻² 之间跳动，应变 0.300 时为 9.6×10¹⁷ m⁻²，HCP 比例稳定在 0.15 到 0.33。跳动来自盒子太小：边长约 36 Å 时，一条贯穿盒子的位错线对应约 8×10¹⁶ m⁻²，位错在周期边界间形核和湮灭会让单帧数值成倍变化。要在论文中比较位错密度，需要更大的盒子，或对一段应变区间取平均。

# plot_ss.py：读取 fix print 输出，画应力-应变曲线，拟合初始模量并找第一次屈服
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

d = np.loadtxt("stress_strain.txt", comments="#")
strain, sxx = d[:, 0], d[:, 1]

# 20 个点滑动平均，压掉热噪声（每点间隔 100 步）；valid 模式避免两端补零造成的假下降
w = 20
sxx_s = np.convolve(sxx, np.ones(w) / w, mode="valid")
strain_s = strain[w // 2 : w // 2 + len(sxx_s)]

lin = strain <= 0.01                     # 只用 0–1% 应变拟合
E = np.polyfit(strain[lin], sxx[lin], 1)[0]

# 第一次屈服：应力超过 1 GPa 后第一次从历史最大值回落 20% 之前的最大值
# 小体系卸载后会再次加载，全局最大值可能出现在后面，不能直接用 argmax
run_max = np.maximum.accumulate(sxx_s)
drop = np.nonzero((run_max > 1.0) & (run_max - sxx_s > 0.2 * run_max))[0]
j = drop[0] if drop.size else len(sxx_s)
i = np.argmax(sxx_s[:j])
print(f"E(0-1%) = {E:.1f} GPa")
print(f"first yield = {sxx_s[i]:.2f} GPa at strain {strain_s[i]:.3f}")

plt.plot(strain, sxx, lw=0.5, alpha=0.4, label="raw")
plt.plot(strain_s, sxx_s, lw=1.5, label=f"moving avg ({w})")
plt.xlabel("Engineering strain")
plt.ylabel("Stress σxx (GPa)")
plt.legend()
plt.tight_layout()
plt.savefig("stress_strain.png", dpi=200)

# 也可以直接解析 log.lammps（需要安装 LAMMPS Python 模块，conda 版自带）
# from lammps.formats import LogFile
# run = LogFile("log.lammps").runs[-1]   # 键名就是 thermo 列名，如 "v_strain"、"v_sxx"

# dxa_density.py：逐帧做 CNA 与 DXA，输出位错密度和 FCC/HCP 原子比例
import sys
from ovito.io import import_file
from ovito.modifiers import DislocationAnalysisModifier

pipeline = import_file(sys.argv[1] if len(sys.argv) > 1 else "tensile.lammpstrj")
dxa = DislocationAnalysisModifier(input_crystal_structure=DislocationAnalysisModifier.Lattice.FCC)
pipeline.modifiers.append(dxa)

print("# step  rho_total(1/m^2)  rho_shockley(1/m^2)  FCC_frac  HCP_frac")
for frame in range(pipeline.source.num_frames):
    data = pipeline.compute(frame)
    a = data.attributes
    vol = a["DislocationAnalysis.cell_volume"]                 # Å^3
    rho = a["DislocationAnalysis.total_line_length"] / vol * 1e20   # Å/Å^3 -> m^-2
    rho_s = a.get("DislocationAnalysis.length.1/6<112>", 0.0) / vol * 1e20
    n = data.particles.count
    fcc = a["DislocationAnalysis.counts.FCC"] / n
    hcp = a["DislocationAnalysis.counts.HCP"] / n
    print(f"{a['Timestep']:8d}  {rho:.3e}  {rho_s:.3e}  {fcc:.3f}  {hcp:.3f}")

续跑

长时间运行与 restart 续跑

restart 文件保存原子、盒子、units、atom_style 等。不保存的内容必须在续跑脚本里重写：所有 fix 和 compute，以及从文件读取参数的多体势（EAM、Tersoff 等）的 pair_coeff。

fix deform 使用和原脚本相同的 fix-ID（这里是 pull）时，会从 restart 文件恢复拉伸开始时的盒子，但只有 run 命令带 start 和 stop 时才使用它。start 填原始拉伸段开始时的步数（reset_timestep 0 之后是 0），stop 填总步数。只写 run 300000 upto 时，deform 以续跑时的盒长为基准重新计算伸长量。本页实测（864 原子、1×10¹⁰ s⁻¹ 的缩短算例，从应变 0.10 续跑到原定终点步）：只写 upto 时终点应变是 0.32；加 start 0 stop 后终点应变是 0.30，终点应力 6.80 GPa，与不中断运行的 6.78 GPa 一致。用本节脚本（只改 L0）从“第一个 in 文件”一节中 4000 原子算例的第 150000 步续跑 1000 步，输出的应变与原运行相同，应力前 7 位有效数字一致。按页面参数从第 150000 步续跑而不加 start/stop，终点应变会是 0.3225。L0 必须填原始拉伸前的长度，不能用续跑时的 lx，否则应变从 0 重新算起。run N upto 让总步数停在 N，中断几次都能对齐到同一终点。

# in.cu_tensile.continue —— 从第 150000 步的 restart 继续拉伸
read_restart    tensile.150000.restart
pair_style      eam/alloy                        # 多体势不存进 restart，必须重写
pair_coeff      * * Cu_mishin1.eam.alloy Cu
neighbor        2.0 bin

variable        T     equal 300.0
variable        erate equal 1.0e9*1.0e-12
variable        L0    equal 72.5                 # 改成原始日志里 variable L0 equal 后面的数值
timestep        0.001
fix             nptyz all npt temp ${T} ${T} 0.1 y 0.0 0.0 1.0 z 0.0 0.0 1.0 drag 1.0
fix             pull all deform 1 x erate ${erate} units box remap x   # fix-ID 与原脚本一致

variable        strain equal (lx-v_L0)/v_L0
variable        sxx equal -pxx/10000
fix             out all print 100 "${strain} ${sxx}" append stress_strain.txt screen no
thermo_style    custom step v_strain temp v_sxx pe press
thermo          1000
restart         50000 tensile.*.restart
run             300000 upto start 0 stop 300000  # 跑到总步数 300000；start/stop 让 deform 沿用原始盒长

报错

新手最常遇到的 4 个报错

Out of range atoms、Too many neighbor bins、NaN 等更多报错的原因和处理见《LAMMPS 常见报错与解决》（/learn/lammps-errors）。

国内环境

国内用户常用做法：镜像、超算 module、离线编译与建模转换

这一节整理计算化学公社、CSDN、知乎和高校超算文档中反复出现的做法，已按 LAMMPS 与相关工具的当前文档或源码核对；标注“本页实测”的数字来自 2026-10-10 的实际查询或运行。

# ~/.condarc：conda-forge 走清华 TUNA 镜像（TUNA 帮助页中只用社区源的写法）
channels:
  - conda-forge
  - nodefaults
show_channel_urls: true
custom_channels:
  conda-forge: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud

conda clean -i
# 构建字符串 = 加速方式 + MPI：cpu、cuda129、cuda130、cuda134 与 nompi、mpich、openmpi 组合
conda create -n lammps -y "lammps=2025.07.22=cpu*openmpi*"          # 新环境总下载约 410 MB
conda create -n lammps-gpu -y "lammps=2025.07.22=cuda129*openmpi*"  # 新环境总下载约 1 GB
conda activate lammps && lmp -h | head -20

# 1) 在能联网的机器上下载第三方源码包（地址与 cmake/Modules/Packages/VORONOI.cmake 中一致）
curl -LO https://download.lammps.org/thirdparty/voro++-0.4.6.tar.gz
# 2) 传到集群后，配置时用 <前缀>_URL 指向本地文件；SHA256 仍按默认值校验，版本必须一致
cd lammps/build
cmake ../cmake -C ../cmake/presets/basic.cmake -D PKG_VORONOI=yes -D DOWNLOAD_VORO=yes \
      -D VORO_URL=file://$HOME/src/voro++-0.4.6.tar.gz
cmake --build . -j 8

# 需要下载的包很多时：联网机器上运行 tools/offline/init_caches.sh 建缓存（默认 ~/.cache/lammps），
# 把缓存目录拷到集群同一路径，再 source tools/offline/use_caches.sh，按它打印的 -D LAMMPS_DOWNLOADS_URL=... 加到 cmake 命令

# VMD Tk Console：Packmol 输出的 pdb 转 LAMMPS data（TopoTools）
package require topotools
mol new system.pdb
# pdb 中 VMD 的 type 取自原子名，C1、C2、H1 会各成一个类型；先按元素合并
[atomselect top {name "C.*"}] set type C
[atomselect top {name "H.*"}] set type H
pbc set {40.0 40.0 40.0}        ;# 没有盒子时 writelammpsdata 拒绝写文件
topo writelammpsdata system.data full
# 类型号按类型名字母序分配：C=1、H=2；Masses 段每行末尾的注释就是类型名

交给 Agent

把这套流程交给 Agent

一句话指令示例：“用 LAMMPS 模拟直径约 4 nm、长约 20 nm 的 [100] 单晶铜纳米线在 300 K、1×10⁹ s⁻¹ 下的单轴拉伸，用 Mishin EAM 势和 GPU 加速，输出应力-应变曲线、杨氏模量、屈服强度和共近邻分析快照，再与直径约 2 nm 的纳米线对比尺寸效应。”

智能体会执行这些步骤：租用 GPU 并运行环境自检；编写建模、弛豫和拉伸的 in 文件；先跑短程测试，估算耗时后再提交正式计算；用 Python 提取应力-应变数据并拟合模量；用共近邻分析标记位错和孪晶；对结果做对抗审阅，确认任务真正完成。交付内容包括 in 文件、后处理脚本、日志、曲线图、结构快照、环境记录和报告。计算结束后释放算力。

以下内容仍需你自己核对：势函数是否适合要研究的变形机制；应力归一化体积的定义；应变速率和尺寸对屈服强度的影响是否已在结论中说明；是否需要用不同随机种子重复计算，以估计统计误差。