LAMMPS 近场动力学 Pair 样式全指南：peri/pmb、peri/lps、peri/ves 与 peri/eps 建模实战

近场动力学（Peridynamics）是经典连续介质力学的非局部推广，特别适合模拟裂纹萌生、扩展与材料碎裂等不连续问题。LAMMPS 通过可选 PERI 包提供了四个近场动力学 pair 样式：peri/pmb（键基原型微弹性脆性模型）、peri/lps（基于状态线性近场动力学固体）、peri/ves（线性粘弹性固体）与 peri/eps（弹塑性固体）。本文以 pair_peri.rst 为骨架，结合 Howto_peri.rst、PERI 包源码 与 examples/peri 示例，系统讲解四种模型的语法、全部 pair_coeff 系数含义、键断裂准则、底层实现原理与工程陷阱，使你能够直接写出可运行、可调参、可后处理的近场动力学输入脚本。

---

一、近场动力学与 LAMMPS PERI 包概述

近场动力学将经典连续介质的运动方程改写为空间积分形式：质点 $\mathbf{x}$ 的运动由邻域 $\mathcal{H}_{\mathbf{x}}$（半径为 horizon $\delta$ 的球形区域）内所有其他质点的力密度之和决定，其运动方程为

$$\rho(\mathbf{x})\,\ddot{\mathbf{u}}(\mathbf{x},t) = \int_{\mathcal{H}_{\mathbf{x}}} \left\{ \underline{\mathbf{T}}[\mathbf{x},t]\langle\mathbf{x}'-\mathbf{x}\rangle - \underline{\mathbf{T}}[\mathbf{x}',t]\langle\mathbf{x}-\mathbf{x}'\rangle \right\} dV_{\mathbf{x}'} + \mathbf{b}(\mathbf{x},t)$$

其中 $\underline{\mathbf{T}}$ 是力矢量态（force vector state），$\mathbf{b}$ 是外体积力密度。horizon $\delta$ 类似于分子动力学中的截断半径。离散化后，每个粒子只与参考构型中距离不超过 horizon 的"成键"粒子相互作用，因此无需经典连续介质理论中的导数与裂纹奇异场假设，断裂自然地通过键的失效产生。

在 LAMMPS 中，PERI 包的核心组件包括（见 src/PERI）：

组件	文件	作用
pair 样式	pair_peri_pmb.cpp、pair_peri_lps.cpp、pair_peri_ves.cpp、pair_peri_eps.cpp	四种本构模型的力与能量计算
基类	pair_peri.cpp、pair_peri.h	公共初始化、影响函数、膨胀量计算
邻接与参考态	fix_peri_neigh.cpp	建立成键伙伴表、参考键长 r0、加权体积 wvolume
原子样式	atom_vec_peri.cpp	管理 vfrac、s0、smin、x0 等逐原子量
诊断计算	compute_damage_atom.cpp、compute_dilatation_atom.cpp	损伤、膨胀量输出

四个样式均属于 PERI 包，需要编译时启用该包（见 Build_package）。peri/pmb 与 peri/lps 还提供 OpenMP 加速变体 peri/pmb/omp、peri/lps/omp（实现位于 src/OPENMP，可通过 -sf omp -pk omp 启用）。

近场动力学的理论基础

近场动力学理论奠基文献为 Silling 2000（键基理论）与 Silling 2007（状态基理论），LAMMPS 中的实现描述见 Parks 2008。在状态基框架下，力矢量态可写为标量力态与变形方向矢量态之积 $\underline{\mathbf{T}} = \underline{t}\,\underline{\mathbf{M}}$，这类材料被称为 ordinary 材料——任意两个相互作用质点间的力沿两点的连线方向作用。peri/lps、peri/ves、peri/eps 均属于状态基模型，而 peri/pmb 是键基模型的经典代表，其力只与键的伸缩量（stretch）线性相关。

---

二、四种模型速览与选型

四个样式分别实现四种本构模型：

• peri/pmb — Prototype Microelastic Brittle（原型微弹性脆性）模型，键基（bond-based）。弹性性质仅由体模量 $K$ 与 horizon $\delta$ 决定；三维下隐含泊松比 $\nu=1/4$（二维为 $1/3$）。键力 $f = c\,s$，其中 $s$ 为键伸缩（见下文），$c$ 为键的弹性常数。
• peri/lps — Linear Peridynamic Solid（线性近场动力学固体）模型，状态基（state-based），是经典各向同性线弹性材料的非局部类比，需要体模量 $K$ 与剪切模量 $G$ 两个参数。
• peri/ves — Linear Viscoelastic Solid（线性粘弹性固体）模型，状态基。在 LPS 弹性响应之上叠加粘弹性松弛项，由松弛参数与时间常数控制。
• peri/eps — Elastic-Plastic Solid（弹塑性固体）模型，状态基。在弹性响应之上叠加理想塑性（屈服应力）行为。

模型核心参数对照表

参数	peri/pmb	peri/lps	peri/ves	peri/eps	含义
c	✓	–	–	–	键弹性常数（能量/距离/体积²）
K	–	✓	✓	✓	体模量（力/面积）
G	–	✓	✓	✓	剪切模量（力/面积）
horizon	✓	✓	✓	✓	相互作用截断距离（距离）
s00	✓	✓	✓	✓	临界伸缩基值（无量纲）
α	✓	✓	✓	✓	临界伸缩压缩修正系数（无量纲）
m_lambdai	–	–	✓	–	粘弹性松弛参数（0~1，无量纲）
m_taubi	–	–	✓	–	粘弹性时间常数（无量纲）
m_yield_stress	–	–	–	✓	屈服应力（力/面积）

---

三、命令语法与 pair_coeff 参数详解

所有四种样式的语法形式一致：

pair_style style

其中 style = peri/pmb 或 peri/lps 或 peri/ves 或 peri/eps（可加 /omp 后缀）。四种样式均不接收额外的全局设置参数。

四种样式的示例（来自 pair_peri.rst 与 examples/peri）：

pair_style peri/pmb
pair_coeff * * 1.6863e22 0.0015001 0.0005 0.25

pair_style peri/lps
pair_coeff * * 14.9e9 14.9e9 0.0015001 0.0005 0.25

pair_style peri/ves
pair_coeff * * 14.9e9 14.9e9 0.0015001 0.0005 0.25 0.5 0.001

pair_style peri/eps
pair_coeff * * 14.9e9 14.9e9 0.0015001 0.0005 0.25 118.43

每种模型通过 pair_coeff 命令为每个原子类型对定义系数，也可以在 read_data / read_restart 读取的数据文件或重启文件中给出。下面逐模型解释。

3.1 peri/pmb：4 个系数

pair_coeff * * c horizon s00 alpha

• c：近场动力学键的有效弹性常数（spring constant），单位量纲为"能量/距离/体积²"。由体模量换算：$c = 18K / (\pi\delta^4)$，因此 c/distance = stiffness/volume²（stiffness 量纲为能量/距离²，volume 为距离³）。示例中 $c = 18\times14.9\times10^9/(\pi(1.5\times10^{-3})^4) \approx 1.6863\times10^{22}$。
• horizon：键的截断距离，决定初始成键范围。
• s00 与 α：键断裂准则参数（见第四节）。

从源码 pair_peri_pmb.cpp 可见，coeff() 严格要求 6 个参数（两个类型范围 + 4 个系数），解析后存入 kspring[i][j]、cut[i][j]、s00[i][j]、alpha[i][j] 数组，并通过 init_one() 对称复制到 [j][i]。

3.2 peri/lps：5 个系数

pair_coeff * * K G horizon s00 alpha

• K：体模量（bulk modulus），量纲为力/面积（Pa）。
• G：剪切模量（shear modulus），量纲为力/面积。
• horizon、s00、α：同 PMB。

LPS 的力标量态为（见 Howto_peri.rst）：

$$\underline{t} = \frac{3K\theta}{m}\,\underline{\omega}\,\underline{x} + \alpha\,\underline{\omega}\,\underline{e}^\mathrm{d}, \qquad \alpha = \frac{15G}{m}$$

其中 $\theta$ 是膨胀量（dilatation），$m$ 是加权体积（weighted volume），$\underline{\omega}$ 是影响函数，$\underline{e}^\mathrm{d}$ 是伸缩标量态的偏量部分。实现中通过 comm_forward = 1 在相邻进程间前向通信膨胀量 $\theta$（见 pair_peri_lps.cpp）。

3.3 peri/ves：7 个系数

pair_coeff * * K G horizon s00 alpha m_lambdai m_taubi

前 5 个参数同 LPS。新增：

• m_lambdai：粘弹性松弛参数，取值在 0 到 1 之间。当 m_lambdai 非常小时，粘弹性模型的响应非常接近线性弹性模型。
• m_taubi：粘弹性时间常数。

源码 pair_peri_ves.cpp 中，粘弹性项按指数衰减实现：

double c1 = taui / timestepsize;   // timestepsize = update->dt
decay = exp(-1.0 / c1);
betai = 1.0 - c1 * (1.0 - decay);

粘弹性回伸量（back extension）$e_{db}^{n+1} = e_d^n(1-decay) + e_{db}^n\cdot decay + betai\cdot\Delta e_d$，总键力为弹性部分与粘弹性部分之和。每步还需将当前偏量伸缩与回伸量存入 fix_peri_neigh 管理的伙伴数组中，因此 VES 比 LPS 多两个逐键状态量。

3.4 peri/eps：6 个系数

pair_coeff * * K G horizon s00 alpha m_yield_stress

前 5 个参数同 LPS。新增：

• m_yield_stress：材料屈服应力，量纲为力/面积。

源码 pair_peri_eps.cpp 中，先对每个粒子计算偏量力态范数 tdnorm 与逐点屈服值 pointwiseYieldvalue = 25.0/(8π) / horizon⁵ · yieldStress²，再判断屈服面函数 fsurf = tdnorm²/2 − pointwiseYieldvalue 是否大于 0；若超过屈服面则进入塑性步，将试探偏量力态投影回屈服面（径向返回式修正），并更新塑性乘子 lambdaValue。

重要版本说明（4Jul2026 起）：早期版本的 peri/eps 塑性模型与文献 Mitchell 2011（SAND2011-3166）相比存在错误，影响塑性偏量伸缩的演化并导致屈服面严重过冲。这些问题已修正，但当前实现仍没有径向返回规则来确保塑性偏量伸缩不离开屈服面，长时间模拟中可能出现偏离屈服面的漂移。这一改进方向记录在 issue #5064 中（源码 pair_peri_eps.cpp 注释亦有说明）。使用 peri/eps 做长时程塑性模拟时应留意该局限。

3.5 系数的单位约定

文档强调，c 的单位使 c/distance = stiffness/volume²。数值算例中（Howto_peri.rst）采用 SI 单位：$K = 14.9$ GPa、horizon $=1.5$ mm、lattice $=0.5$ mm。若使用其他单位体系，请按上述量纲关系换算。

---

四、键断裂准则（Bond Breaking Criterion）

4.1 临界伸缩公式

四种样式共用同一套键断裂准则。键伸缩定义为

$$s = \frac{r - r_0}{r_0}$$

其中 $r$ 是当前两粒子间距离，$r_0$ 是参考构型中的初始键长（键伸缩等价于一维应变：平衡键长时 $s=0$，两倍键长时 $s=1$）。当键伸缩超过临界伸缩时，键不可逆地断裂。

根据 Parks 2008（eq. 9），粒子的临界伸缩为

$$s_0 = s_{00} - \alpha\, s_{\min}$$

其中 $s_{\min}$ 是该粒子所有键中最小（最受压）的伸缩量。一条键断裂的条件是 $s > \min(s_{0,i},\, s_{0,j})$，即从两端粒子对称地判断——两个粒子必须对它们的共同键采用相同的判定。

4.2 逐键评估（4Jul2026 起的关键变更）

s00 与 α 可以为每种原子类型对单独赋值，例如给两种材料界面处的键设定更小的 s00，使裂纹优先在界面萌生。临界伸缩逐键使用该键自己的 s00 与 α，结合每个端点的几何 $s_{\min}$ 进行求值。

早期版本将临界伸缩存储为单个逐粒子量（取粒子所有键中 $s_{00}-\alpha s$ 的最大值），这只有在所有类型对的 s00、α 完全相同时才等价于上述准则；一旦系数随类型对变化，键会错误地继承周围体相键的临界伸缩，导致弱化界面无法断裂。4Jul2026 起该准则改为逐键评估。

源码中的实现（以 pair_peri_pmb.cpp 为例）：

if (stretch > s00[itype][jtype] - alpha[itype][jtype]*MAX(smin[i],smin[j]))
  partner[i][jj] = 0;   // 键不可逆断裂

其中 stretch = dr / r0[i][jj]，smin[i] 与 smin[j] 是上一个时间步存储的逐粒子最小伸缩。为效率考虑，实现不在一趟内先算 $s_{\min}$ 再判断，而是沿用上一步的 $s_{\min}$（详见 Howto_peri.rst 的 Damage 一节）；数学上 $\min(s_{0,i},s_{0,j}) = s_{00}-\alpha\,\max(s_{\min,i},s_{\min,j})$，与逐键判据一致。

compute property/atom 的 s0 逐原子属性报告的就是每个粒子的临界伸缩。早期文档曾错误地把 s0 描述为粒子参与任何键的最大伸缩，现已更正（可对照 atom_vec_peri.cpp 中 s0、smin 的初始化）。

4.3 首个时间步无断键

技术细节上，键直到第一个时间步结束才建立完成，因此不应在模拟的前两步尝试断键。atom_vec_peri.cpp 中 smin 初始化为 -DBL_MAX（使隐含临界伸缩为 $+\infty$），pair_peri_pmb.cpp 在每步末尾对无幸存键的原子保持 -DBL_MAX 哨兵值，保证它不会通过 max() 触发邻居键的断裂。在数值算例中，弹丸被放置于几步之后才接触靶板，正是为了避开这一约束。

---

五、从零搭建近场动力学模拟

5.1 最小可运行输入（Quick Start）

以下最小示例来自 Howto_peri.rst：

units         si
boundary      s s s
lattice       sc 0.0005
atom_style    peri
atom_modify   map array
neighbor      0.0010 bin
region        target cylinder y 0.0 0.0 0.0050 -0.0050 0.0 units box
create_box    1 target
create_atoms  1 region target

pair_style    peri/pmb
pair_coeff    * * 1.6863e22 0.0015001 0.0005 0.25
set           group all density 2200
set           group all volume 1.25e-10
velocity      all set 0.0 0.0 0.0 sum no units box
fix           1 all nve
compute       1 all damage/atom
timestep      1.0e-7

逐行要点：

• units si：近场动力学模拟通常使用 SI 单位。
• atom_style peri：必须使用 PERI 原子样式，它管理 rmass、vfrac、s0、smin、x0 五个逐原子量（见 atom_vec_peri.cpp）。
• atom_modify map array：必须建立原子映射，用于按 tag 检索成键伙伴。
• lattice sc 0.0005：粒子必须生成在简单立方（simple cubic）格点上；源码 pair_peri.cpp 还会强制校验 x/y/z 三个方向的晶格常数完全一致，否则报错。
• neighbor 0.0010 bin：skin 距离需满足"horizon + skin > 最大可能成键距离"（见第八节陷阱）。这里取 0.001 m。
• set ... density 2200 / set ... volume 1.25e-10：必须为每个粒子显式定义质量密度与体积分额。简单立方格点上粒子体积应等于晶格常数的立方，即 $V_i = \Delta x^3 = (0.0005)^3 = 1.25\times10^{-10}$。
• fix 1 all nve：速度-Verlet 积分器（等价于时间中心差分，数值上更稳定）。
• compute 1 all damage/atom：逐原子损伤诊断（见第七节）。

5.2 刚性球冲击脆性圆盘：完整数值算例

这是复现 Silling 2005 论文第 6 节实验的经典算例（完整注释脚本见 Howto_peri.rst 与 examples/peri/in.peri.pmb）：

问题设定：直径 0.01 m、速度 100 m/s 的刚性球垂直冲击均质脆性圆盘。靶材密度 2200 kg/m³，PMB 模型 $K = 14.9$ GPa，$s_{00}=0.0005$、$\alpha=0.25$。三维简单立方格点晶格常数 0.0005 m、horizon 0.0015 m（3 倍晶格常数）。靶为直径 0.074 m、厚度 0.0025 m 的圆柱，共 103,110 个粒子，每个粒子体积 $V_i=1.25\times10^{-10}$ m³。CFL 分析表明 $1.0\times10^{-7}$ s 的时间步长是安全的。

# 3D Peridynamic simulation with projectile
units           si
dimension       3
boundary        s s s
atom_style      peri
atom_modify     map array
neighbor        0.0010 bin
lattice         sc 0.0005
# Create desired target
region          target cylinder y 0.0 0.0 0.037 -0.0025 0.0 units box
# Make 1 atom type
create_box      1 target
# Create the atoms in the simulation region
create_atoms    1 region target
pair_style      peri/pmb
#               <type1> <type2>    <c>    <horizon>  <s00> <alpha>
pair_coeff         *       *    1.6863e22 0.0015001 0.0005  0.25
# Set mass density
set             group all density 2200
# volume = lattice constant^3
set             group all volume 1.25e-10
# Zero out velocities of particles
velocity        all set 0.0 0.0 0.0 sum no units box
# Use velocity-Verlet time integrator
fix             F1 all nve
# Construct spherical indenter to shatter target
variable        y0 equal 0.00510
variable        vy equal -100
variable        y equal "v_y0 + step*dt*v_vy"
fix             F2 all indent 1e17 sphere 0.0000 v_y 0.0000 0.0050 units box
# Compute damage for each particle
compute         C1 all damage/atom
timestep        1.0e-7
thermo          200
dump            D1 all custom 100 dump.peri id type x y z c_C1
run             2000

要点说明：

• 弹丸不是 Silling 2005 中的原始弹丸，而是施加 $F(r) = -k_s(r-R)^2$ 的球面压头（fix indent），$R = 0.05$ m、$k_s = 1.0\times10^{17}$；力始终排斥，$r>R$ 时 $F=0$。
• 换用 LPS 模型只需把 pair_style 换成 peri/lps、pair_coeff 换成 14.9e9 14.9e9 0.0015001 0.0005 0.25；examples/peri 目录提供了 pmb/lps/ves/eps 四个完整可运行脚本（in.peri.pmb、in.peri.lps、in.peri.ves、in.peri.eps），VES 与 EPS 的 pair_coeff 分别追加 0.5 0.001 与 10.0e8。
• 浮点精度提示：IEEE 双精度下，未伸缩键的 fl(s) 可能为约 $2.22\times10^{-16}$ 的机器精度量级残差，导致 $c\,s\,V_i \approx 4.68\times10^{-4}$ 的虚假力。源码中通过 if (fabs(dr) < NEAR_ZERO) dr = 0.0;（NEAR_ZERO = 2.2204e-16，定义于 pair_peri.h）消除该误差。

---

六、源码实现原理纵深：力、短程接触与体积修正

6.1 通用结构：固定伙伴表 + 双计数

所有 peri pair 样式共享相同计算骨架（pair_peri.cpp）：init_style() 强制要求 atom_style peri、原子映射与简单立方晶格，并自动添加内部 fix PERI_NEIGH all PERI_NEIGH。fix_peri_neigh.cpp 在第一次运行时扫描完整邻域表，将 horizon 内的粒子对写入 partner 伙伴表，记录参考键长 r0<a href="https://link.gitcode.com/i/c49d01210d6b366cc6ab52f082695a05" target="_blank">i][jj]、中间体积 vinter 与加权体积 wvolume，并在日志中打印键统计（"Peridynamic bonds: total # of bonds / bonds per atom"）。伙伴表在进程迁移、重启文件中会被压缩打包（[fix_peri_neigh.cpp），pair_peri_pmb.cpp 的伙伴循环中"partner 列表含每条键两次（I-J 出现两遍）"，因此能量与应力统计使用 1/2 因子并以 (i,i) 记账。

周期边界下若模拟盒小于 2 倍 horizon，同一个粒子可能在邻域表中出现两次，fix_peri_neigh 会报错"Duplicate particle in PeriDynamic bond - simulation box is too small"（fix_peri_neigh.cpp）。

6.2 影响函数与加权体积

状态基模型（LPS/VES/EPS）依赖影响函数 $\underline{\omega}$。LAMMPS 默认实现为

$$\underline{\omega}\langle\|\boldsymbol{\xi}\|\rangle = 1/\|\boldsymbol{\xi}\|$$

见 pair_peri.h 的 influence_function()。用户可通过修改 pair_peri_lps.cpp 中的同名方法自定义影响函数；代码同时支持球面与非球面影响函数（各向同性/各向异性材料）。

加权体积

$$m[\mathbf{x}] = \int_{\mathcal{H}_\mathbf{x}} \underline{\omega}\langle\boldsymbol{\xi}\rangle\, \underline{x}\langle\boldsymbol{\xi}\rangle^2\,dV_{\boldsymbol{\xi}}$$

在 fix_peri_neigh.cpp 中按参考构型逐键累加（PMB 使用 $\omega=1$），并在每次重排邻域时前向通信给 ghost 原子。注意加权体积 $m$ 与时间无关，即使键断裂也不改变——它始终基于参考构型的键族计算。

6.3 膨胀量的通信与能量

LPS/VES/EPS 每步先调用 compute_dilatation()（pair_peri.cpp）累加 $\theta$，再通过 comm->forward_comm(this) 将膨胀量广播到 ghost 原子（comm_forward = 1）。体积相关能量项为 $0.5\,K\,\theta^2$，在 pair_peri_lps.cpp 中分别累加到全局能与逐原子能。

6.4 短程接触力

为避免粒子在键断裂或强烈压缩时互相穿透，四种样式都引入了纯排斥的短程力（pair_peri_pmb.cpp）：

• 短程相互作用距离 $d_{ij} = \min\{0.9\,\|\mathbf{x}_j-\mathbf{x}_i\|,\; 1.35\,l_c\}$（$l_c$ 为晶格常数）；
• 接触刚度 $c_S = 15\,c = 15\times 18K/(\pi\delta^4)$（状态基模型用 $c_S = 15\cdot 18K/(\pi\,\delta^4)$，见 pair_peri_lps.cpp）；
• 仅当 $r < d_{ij}$ 时产生排斥力，$dr = r - d_{ij} < 0$，力方向沿粒子连线。

短程力只在粒子受压时出现，这保证了碎片间的接触响应。

6.5 节点体积修正（horizon 边界缩放）

位于 horizon 边缘附近的粒子，其体积并未被完全包含在 horizon 内，若直接用 $V_j$ 会高估对方看到的作用体积。实现采用线性无量纲缩放函数 $\nu(\mathbf{x}_p - \mathbf{x}_i)$（pair_peri_pmb.cpp）：

• 若 $\delta - r_i \le \|\mathbf{x}_p-\mathbf{x}_i\| \le \delta$：$\nu$ 从 1 线性降至 0.5（在 $\|\cdot\|=\delta$ 时 $\nu=0.5$）；
• 若 $\|\mathbf{x}_p-\mathbf{x}_i\| \le \delta - r_i$：$\nu=1$；
• 否则为 0。

其中 $r_i$ 取晶格常数的一半。该修正同时作用于键力、加权体积与膨胀量计算。

---

七、损伤诊断与可视化

7.1 损伤计算

损伤（damage）定义为某粒子周围键失效的比例（Silling 2005）：

$$\varphi(\mathbf{x}, t) = 1 - \frac{\int_{\mathcal{H}_\mathbf{x}} \mu(t,\boldsymbol{\eta},\boldsymbol{\xi})\,dV_{\mathbf{x}'}} {\int_{\mathcal{H}_\mathbf{x}} dV_{\mathbf{x}'}}$$

其中 $\mu$ 是历史相关布尔函数（键存活则 1）。LAMMPS 中由 compute damage/atom 实现（compute_damage_atom.cpp）：遍历粒子的伙伴表，累加存活伙伴的体积 vfrac<a href="https://link.gitcode.com/i/5421ad431fca41a19e61b98712d8209c" target="_blank">j]，damage = 1 - 存活体积/初始中间体积 vinter。数值在 0（无损伤）到 1（所有键断裂）之间。注意该计算依赖内部 fix PERI_NEIGH，若在非近场动力学势下使用会报错（[compute_damage_atom.cpp）。

compute dilatation/atom（compute_dilatation_atom.cpp）则通过 pair 的 extract("theta") 接口输出每个粒子的膨胀量，但不能与 peri/pmb 一起使用（PMB 无膨胀量概念，compute_dilatation_atom.cpp 明确报错）。

7.2 用 dump 生成图像与动画

可以直接用 dump image / dump movie 将损伤着色的粒子快照输出为图片序列：

dump            D2 all image 100 dump.peri.*.png c_C1 type box no 0.0 view 30 60 zoom 1.5 up 0 0 -1  ssao yes 4539 0.6
dump_modify     D2 pad 5 adiam * 0.001 amap 0.0 1.0 ca 0.1 3 min blue 0.5 yellow max red

若需要更精细的交互式可视化，可先用 dump custom 写出 id type x y z c_C1 轨迹（数值算例中的 dump D1 ... dump.peri），再导入 OVITO：加载 dump.peri → 将粒子显示半径设为 0.0005 → 添加 "Color coding" 修饰（颜色梯度 Jet，Start=0、End=1），必要时叠加 "Ambient occlusion" 提升立体感。相关运行截屏与示例图位于 doc/src/JPG（如 dump.peri.300.png、dump.peri.600.png、dump.peri.2000.png、ovito-peri-snap.png）。

7.3 关于解对称性

采用"完美"格点与正对中心的完美球形弹丸时，计算结果呈现明显的对称性。若需打破对称，可用 displace_atoms 命令对格点位置做微小扰动。

---

八、BPM 框架下的替代实现（新选择）

pair_peri.rst 特别提醒：BPM（bonded particle model）包提供了上述四种模型的另一种实现——bond_style bpm/peri 加接触 pair_style bpm/peri。它将近场动力学键变成真正的 LAMMPS 键，使用标准 atom_style bond（无需专用原子样式），可以用 create_bonds / delete_bonds 或数据文件控制实体域与预置裂纹、空洞等缺陷，且显著快于本文档所述的 peri pair 样式。两者遵循相同的本构方程（解析弹性响应一致）；据 Howto_peri.rst 的基准，在周期 16×16×16 简单立方块（4096 节点）上，BPM 实现每步力计算的耗时约为 legacy peri/* 的 1/3.6（pmb）、1/5.4（lps）、1/4.4（ves）与 1/5.9（eps）（单 MPI rank 下，4 rank 约低 4 倍），主要原因是 legacy 样式会双计数遍历伙伴表，而 BPM 每条键只访问一次。详细的 PERI→BPM 迁移对照表（atom_style peri→atom_style bond、set volume→set d_vfrac、compute damage/atom→compute property/atom d_damage 等）见 Howto_peri.rst。如何取舍可参阅 Howto_peri.rst 的逐模型耗时对比。

---

九、混合、修饰选项、重启与 rRESPA

pair_peri.rst 对四种样式统一声明了以下行为：

• 不支持 mixing：所有 I,J 类型对的系数都必须显式指定，不能依赖自动混合规则。
• 不支持 pair_modify shift 选项。
• pair_modify table 与 tail 选项不适用。
• 支持二进制重启文件：四种样式都会把自己的系数写入二进制重启文件（write_restart / read_restart，见各源文件末尾），因此读取重启文件的输入脚本无需重复 pair_style / pair_coeff 命令。内部 fix PERI_NEIGH 同样支持全局与逐原子重启（fix_peri_neigh.cpp），重启后可继续原本的键状态。
• rRESPA：只能在 run_style respa 的 pair 关键字下使用，不支持 inner、middle、outer 关键字。
• 默认值：无（所有系数必须显式给出）。

---

十、工程陷阱与调试要点

10.1 skin 距离过小导致崩溃

neighbor 命令的 skin 决定邻域表保留的粒子范围（截断 = horizon + skin）。若 skin 太小，成键粒子在键断裂前可能跑出邻域表，导致程序意外崩溃。经验法则：skin 应使 horizon + skin 大于两个成键粒子的最大可能间距。例如增大 $s_{00}$ 允许更大伸缩时，可能需要同步增大 skin（Howto_peri.rst 的 Pitfalls 一节）。

10.2 horizon 的定义精度

建议将 horizon 设为"略大于整数倍晶格常数"的值。例如晶格常数 0.0005、希望 horizon 为 3 倍晶格常数时，应设 $\delta = 0.0015001$ 而非 0.0015——否则浮点误差可能导致部分相距 3 倍晶格常数的粒子对不成键。

10.3 "Lost" 粒子

若粒子漂移出模拟盒即成为 lost 粒子。默认行为是报错退出；可通过 thermo_modify lost ignore 忽略或 lost warn 告警。近场动力学模块自身不检查 lost 粒子——如果携带未断键的粒子丢失，其键会被剩余粒子标记为断裂。碎片化模拟中建议用固定（非 shrink-wrap）边界并配合 thermo_modify lost ignore（BPM 实现中则是显式、可选的行为，见 Howto_peri.rst）。

10.4 并行扩展性

LAMMPS 按空间分解并行，各子域大致等体积分配。近场动力学模拟中常出现部分子域粒子密集、部分稀疏的负载不均衡，影响并行扩展。这是近场动力学类问题的固有特征，规划大规模并行计算时应预估。

10.5 测试验证：单元测试如何确认行为

PERI 包配有独立的单元测试 unittest/force-styles/test_pair_peri.cpp，自包含地验证物理不变量：

• pmb_interface_fracture_984 / lps_interface_fracture_984：构造双材料 bar，界面处 s00 取小值，验证弱界面键优先断裂（这正是 4Jul2026 逐键判据修复对应的回归测试）；
• no_break_below_threshold：低于临界伸缩时无键断裂（验证 smin 首步初始化不产生虚假断裂）；
• s0_equals_critical_stretch：验证诊断量 $s_0 = s_{00} - \alpha s_{\min}$；
• momentum_conservation：fix nve 无外力下总动量守恒；
• omp_consistency：peri/pmb/omp、peri/lps/omp 加速变体与基础样式逐原子力、损伤一致。

这些测试是修改或扩展近场动力学模块后的回归保障。开发新模型时，Howto_peri.rst 建议：键基材料以 peri/pmb 为模板，状态基材料以 peri/lps 为模板（涉及膨胀量通信与加权体积的部分可直接复用基类）。

---

参考文献

• (Silling 2000) Silling, J Mech Phys Solids, 48, 175-209 (2000) — 键基近场动力学理论。
• (Silling 2007) Silling, Epton, Weckner, Xu, Askari, J Elasticity, 88, 151-184 (2007) — 状态基理论。
• (Parks 2008) Parks, Lehoucq, Plimpton, Silling, Comp Phys Comm, 179(11), 777-783 (2008) — LAMMPS 中的近场动力学实现（临界伸缩准则 eq. 9）。
• (Mitchell 2011) Mitchell, A non-local, ordinary-state-based viscoelasticity model for peridynamics, Sandia National Lab Report 8064 (2011) — VES 模型。
• (Mitchell 2011a) Mitchell, A Nonlocal, Ordinary, State-Based Plasticity Model for Peridynamics, Sandia National Lab Report 3166 (2011) — EPS 模型。
• (Silling 2005) Silling, Askari, Computers and Structures, 83, 1526-1535 (2005) — 数值算例来源。
• (Seleson 2010) Seleson, Parks, Int J Mult Comp Eng, 9(6), 689-706 (2011) — 影响函数与 PMB/LPS 关系。

关于 VES/EPS 的更多实现细节，LAMMPS 发行版还附带 doc/PDF/PDLammps_VES.pdf 与 doc/PDF/PDLammps_EPS.pdf 两份补充文档。发现近场动力学模块缺陷时，可在 GitHub Issue 中按最小复现原则提交（含最少原子数、最少进程数与最简输入脚本），并附问题描述与代码定位线索。