企业官网建设流程全解析

• 科研
• 科学计算
• 高性能计算

本文以 LAMMPS 官方文档 pair_meam_sw_spline.rst 为核心骨架，结合 src/MANYBODY/pair_meam_sw_spline.cpp 源码实现、potentials 目录下的势文件以及 examples/PACKAGES/meam_sw_spline 示例脚本，系统讲解 LAMMPS 中meam/sw/spline对势样式的物理模型、势文件格式、pair_coeff 配置方法、底层实现原理与运行限制。读完本文，你将能够独立完成 Si、Ti 等金属体系的 MEAM+SW 势配置、内聚能验证与能量守恒测试，并理解其在 MANYBODY 包中的调用机制。

一、MEAM+SW：把修正嵌入原子法与 Stillinger-Weber 三项结合

meam/sw/spline是 LAMMPS 中用于金属体系的一种多体对势样式（pair style），它在**修正嵌入原子法（Modified Embedded Atom Method, MEAM）**的基础上（参考 Lenosky 等人 2000 年的工作）额外叠加了一个Stillinger-Weber（SW）三项角项（参考 Stillinger & Weber 1985 年的经典工作），因此被称为 MEAM+SW。该势格式最早由俄亥俄州立大学的 Jeremy Nicklas、Michael Fellinger 和 Hyoungki Park 提出，用于开发 bcc 过渡金属势（这一点在 pair_meam_sw_spline.cpp 文件末尾的版权注释中也有明确记载）。

1.1 总能量分解

MEAM+SW 势的总能量 E 由两部分构成：

E = E_MEAM + E_SW E_MEAM = Σ_IJ φ(r_IJ) + Σ_I U(ρ_I) E_SW = Σ_I Σ_JK F(r_IJ) F(r_IK) G(cos(θ_JIK)) ρ_I = Σ_J ρ(r_IJ) + Σ_JK f(r_IJ) f(r_IK) g(cos(θ_JIK))

其中：

• ρ_I 是原子 I 处的（电子）密度；
• θ_JIK 是以原子 I 为中心、J 与 K 之间的夹角；
• 七个函数φ、F、G、U、ρ、f、g全部用**三次样条（cubic spline）**表示。

从物理图像上看：E_MEAM 部分继承了经典 EAM 的"对势 φ + 嵌入能 U(ρ)"结构，但密度 ρ 中除了二体贡献 Σ_J ρ(r_IJ) 外，还加入了三项角贡献 Σ_JK f(r_IJ) f(r_IK) g(cos θ)，这使势具备方向性成键能力；E_SW 部分则完全由三项角贡献 F(r_IJ)F(r_IK)G(cos θ) 组成，相当于显式加入 Stillinger-Weber 式的键角项。二者叠加后，可以同时描述金属键合与方向性共价键特征。

1.2 与 meam、meam/spline 的关系

在 LAMMPS 中，meam/sw/spline属于 MANYBODY 包中 MEAM 家族的一员，与 pair_style meam、pair_style meam/spline 相关但实现思路不同：

• pair_style meam：经典 MEAM，势函数以解析形式描述；
• pair_style meam/spline：样条化的 MEAM（无 SW 三项叠加项，且支持多组分）；
• pair_style meam/sw/spline：样条化 MEAM 与 SW 角项的叠加，目前仅支持单元素。

在 doc/src/pair_style.rst 中，meam/sw/spline的定位是 "Splined version of MEAM with a Stillinger-Weber term"（带 Stillinger-Weber 项的样条化 MEAM）。

二、语法与最小可用配置

2.1 pair_style 语法

meam/sw/spline的 pair_style 命令不接受任何参数：

pair_style meam/sw/spline

这一点在源码中得到了确认：pair_meam_sw_spline.cpp 的settings()函数规定：

void PairMEAMSWSpline::settings(int narg, char **/*arg*/) { if (narg != 0) error->all(FLERR,"Illegal pair_style command"); }

即任何多余参数都会直接报错Illegal pair_style command。

2.2 官方文档示例

pair_style meam/sw/spline pair_coeff * * Ti.meam.sw.spline Ti # 3 个原子类型全部映射为 Ti 元素： pair_coeff * * Ti.meam.sw.spline Ti Ti Ti

第一条 pair_coeff 适用于单类型体系；第二条示例假设模拟中有 3 个 LAMMPS 原子类型且全部用 Ti 势处理，则文件名后跟 3 个 Ti 元素名完成映射。

2.3 一个可运行的完整示例（内聚能验证）

仓库 examples/PACKAGES/meam_sw_spline/Si/in.fcc 给出了最简的 Si fcc 结构验证脚本，完整命令如下：

units metal boundary p p p atom_style atomic lattice fcc 4.147 region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box pair_style meam/sw/spline pair_coeff * * Si.b.meam.sw.spline Si mass * 28.085 variable cohesive_energy equal pe/atoms variable atmVol equal vol/atoms variable aLatt equal (4*vol/atoms)^0.3333333333 run 0 print "Calculated cohesive energy: ${cohesive_energy} eV/atom" print "Reference cohesive energy: -4.288 eV/atom" print "Atomic volume ${atmVol} A^3" print "Lattice constant ${aLatt} A" print "Reference lattice constant 4.147 A"

Ti 的 hcp（α 相）验证脚本见 examples/PACKAGES/meam_sw_spline/Ti/in.alpha，其中使用lattice custom构建 hcp 晶格，参考内聚能为-4.831 eV/atom，质量设置为mass * 47.90。

2.4 三个必须注意的配置细节

1. 截断半径在势文件里，不在命令里：与其他势不同，MEAM+SW 的截断不能通过 pair_style 或 pair_coeff 命令设置，而是由势文件本身定义。从源码看，read_file() 在读取完势文件后，取 φ、ρ、f、F 四个函数样条截断的最大值作为全局cutoff，并写入cutsq[i][j]供邻表构建使用。
2. 必须用 mass 命令指定原子质量：与 EAM 从势文件读取质量不同，.meam.sw.spline势文件不包含质量信息，因此必须显式执行mass命令。示例中为mass * 28.085（Si）或mass * 47.90（Ti）。
3. 参数基于 metal 单位：所有.meam.sw.spline势文件均按 LAMMPS metal 单位（能量 eV、距离 Å）标定，因此输入脚本中应使用units metal。势文件首行注释中的UNITS: metal字段也明确标注了这一点。

三、势文件格式深入解析

.meam.sw.spline势文件分布在仓库的 potentials 目录（当前仓库提供 Si 与 Ti 两种元素），同时 examples/PACKAGES/meam_sw_spline 下也各有一份副本。以 potentials/Si.b.meam.sw.spline 为例：

# DATE: 2012-10-26 UNITS: metal CONTRIBUTOR: Robert Rudd, robert.rudd@llnl.gov CITATION: Balamane et al, Phys Rev B 46, 2250 (1992) COMMENT: ... 37 -5.511362311759948e+02 0.000000000000000000e+00 1 0 1 0 9.2969999999999997E-01 1.1344841447915871E+02 2.9969523122045521E+03 ...

文件的组织方式如下：

在 potentials/Ti.meam.sw.spline 中可以看到同样 7 段的布局（节点数依次为 13、2、2、11、4、10、8）。每行节点数据中的第三列y2是预计算的二阶导数，源码解析时会将其忽略（见 parse() 中// double y2 = values.next_double(); ignored的注释），二阶导数由程序在 prepareSpline() 中基于自然样条边界条件自行求解三对角方程得出。

从源码结构看，样条函数在内部由SplineFunction类管理（定义于 pair_meam_sw_spline.h），每个样条持有节点位置 X、节点值 Y、二阶导数 Y2、一阶差分 Ydelta 以及两端导数 deriv0/derivN。特别值得说明的是宏SPLINE_MEAMSW_SUPPORT_NON_GRID_SPLINES：

• 置 0（当前仓库默认）：要求样条节点落在均匀网格上，插值时可通过klo = (int)(x * inv_h)直接定位区间，代码注释声称可提速约 50%；
• 置 1：支持非均匀节点，但需二分查找区间，速度较慢。

prepareSpline()会在读取节点时检查fabs(h*i + xmin - X[i]) > 1e-8来判断是否均匀网格；若节点非均匀而宏为 0，会抛出 TokenizerException 报错。也就是说，使用自定义非均匀网格势文件前，需要重新编译时开启该宏。

四、pair_coeff 的元素映射规则

meam/sw/spline仅使用一条pair_coeff 命令，其参数格式为：

pair_coeff * * <势文件名> <N 个元素名>

其中 N 等于 LAMMPS 原子类型数，元素名依次映射到原子类型 1、2、…、N：

1. 前两个参数必须是* *，以覆盖全部 LAMMPS 原子类型；
2. 势文件名指定包含所有所需元素参数的势文件；
3. N 个元素名按顺序将 LAMMPS 原子类型映射到势文件中的元素。

以 3 个原子类型全部为 Ti 为例：

pair_coeff * * Ti.meam.sw.spline Ti Ti Ti

即原子类型 1、2、3 均映射为势文件中的 Ti 元素。若映射值为NULL，则该原子类型不执行映射——这主要用于pair_style hybrid场景：NULL 作为占位符，表示该原子类型将交由其他势处理。例如一个 hybrid 脚本中，类型 1 用 meam/sw/spline，类型 2 用其他对势：

pair_style hybrid meam/sw/spline lj/cut 2.5 pair_coeff * * meam/sw/spline Ti.meam.sw.spline Ti NULL pair_coeff * * lj/cut 1 1 0.1 3.0

关于势文件的搜索路径，可参考 pair_coeff 页面对路径的替代指定方式（如pair_coeff中直接使用路径或通过 LAMMPS_POTENTIALS 环境变量定位）。

源码层面对元素映射的限制在 coeff() 中：

map_element2type(narg-3,arg+3); // for now, only allow single element if (nelements > 1) error->all(FLERR, "Pair meam/sw/spline only supports single element potentials");

即当前实现只允许映射到单一元素，若出现两个不同元素名会直接报错Pair meam/sw/spline only supports single element potentials。这意味着：

• 原子类型可以多于 1（如Ti Ti Ti），但所有类型必须映射到同一个元素；
• 合金/多组分体系目前不支持（官方文档同样明确：currently supports only single-element MEAM+SW potentials. It may be extended for alloy systems in the future）。

五、源码级实现原理：一次 compute() 的三阶段旅程

理解meam/sw/spline的底层执行流程，有助于诊断能量异常或性能问题。从 compute() 的实现看，每次力与能量计算分为三个阶段：

阶段 1：遍历全邻表（full neighbor list），累计密度并计算嵌入能

对每个原子 I，遍历其全部邻居 J：

• 用f.eval(rij)、F.eval(rij)求二体量；
• 对每对已记录的键 K，计算cos θ并累计partial_sum += f_K * g(cosθ)与partial_sum2 += F_K * G(cosθ)（见 pair_meam_sw_spline.cpp），对应公式中的三项密度项；
• 同时累计二体密度rho.eval(rij)；
• 随后用U.eval(rho_value, Uprime_i)求嵌入能 U(ρ) 及其导数 U′(ρ)，并扣除真空零点能zero_atom_energy（其值为U.eval(0.0)，见 read_file()，目的是让真空中孤立原子的能量为 0）。

阶段 2：forward_comm 通信，交换 U′(ρ)

嵌入能导数 U′(ρ_I) 需要通过 MPI 前向通信广播给邻居（pack_forward_comm/unpack_forward_comm），为阶段 3 的二体对力做准备。

阶段 3：遍历半邻表（half neighbor list），计算二体对势与最终力

对每个原子对 I-J：

• fpair = rho_prime(rij) * (Uprime_values[i] + Uprime_values[j]) + phi_prime(rij)，即二体对势 φ 的导数与嵌入能导数贡献之和；
• 除以 r_ij 得到梯度方向的力并累加到 I、J 上（pair_meam_sw_spline.cpp）；
• 三项贡献的力已在阶段 1 的三体循环中通过ev_tally3(i, j, k, ...)计入（pair_meam_sw_spline.cpp）。

也正因这种"全邻表算密度、半邻表算对力"的算法结构，init_style() 中同时请求了 FULL 和 half 两张邻表，并强制要求newton pair on（newton off 时报错Pair style meam/sw/spline requires newton pair on）。

六、验证方法：单元测试与示例基准

6.1 官方回归测试基准

仓库 unittest/force-styles/tests/atomic-pair-meam_sw_spline.yaml 提供了该 pair style 的自动回归测试配置：

lammps_version: 2 Sep 2026 epsilon: 2e-12 prerequisites: | pair meam/sw/spline pre_commands: | atom_modify map array units metal lattice diamond 5.43 region box block 0 2 0 2 0 2 create_box 1 box create_atoms 1 box displace_atoms all random 0.1 0.1 0.1 623426 mass 1 29.0 velocity all create 100 4534624 pair_style: meam/sw/spline pair_coeff: | * * Si.b.meam.sw.spline Si

该测试以 64 个 Si 原子（diamond 5.43 晶格 + 随机位移）为基准，给出了精确到小数点后 14 位的初始化能量（init_vdwl: -110.45990817806636）与运行后能量（run_vdwl: -111.18747862613623），并指定 2e-12 的能量容差。读者可将其作为自建势文件正确性的校准基准：若用同一输入重算得到一致的能量与力，则说明势文件解析与样条插值均正常。

6.2 真空中单原子测试

examples/PACKAGES/meam_sw_spline/Si/in.single_atom 用固定边界把单个 Si 原子放在真空中（boundary f f f，box 范围 -100～100 Å）：

pair_style meam/sw/spline pair_coeff * * Si.b.meam.sw.spline Si mass * 28.085 variable cohesive_energy equal pe/atoms run 0

输出参考内聚能为0 eV/atom。这验证了前述zero_atom_energy = U.eval(0.0)的"零点能扣除"机制：孤立原子的能量恰好被归零。

6.3 NVE 能量守恒测试

examples/PACKAGES/meam_sw_spline/Si/in.energy_conservation.meam.sw 构建 diamond 5.431 的 Si 块体，用fix nve在 0.002 ps 时间步下运行 2000 步，并以thermo_style custom step vol etotal press pxx pyy pxz监控总能量是否守恒——这是检验多体势力场正确性的经典手段（Ti 版本见 examples/PACKAGES/meam_sw_spline/Ti/in.energy_conservation.meam.sw，其中还演示了用lattice custom构建六方晶格的方法）。

七、兼容性与使用限制汇总

7.1 混合（mixing）、shift、table、tail

• meam/sw/spline不支持多种元素类型间的混合（mixing），本身即为单元素势；
• 不支持pair_modify 的shift、table、tail选项。

7.2 重启文件与 rRESPA

• 该势的信息不写入二进制重启文件（因为完全存储在外部势文件中），因此读取 restart 文件后必须在输入脚本中重新指定 pair_style 与 pair_coeff；
• 配合 run_style respa 时，只能通过pair关键字使用，不支持inner、middle、outer关键字。

7.3 运行环境与构建要求

• 必须开启newton on（newton 命令），pair 相互作用要求 Newton 模式开启；
• 该 pair style 仅当 LAMMPS 以MANYBODY 包构建时才可用。CMake 方式下需启用：

cmake -D PKG_MANYBODY=on ../cmake

或传统 make 方式下在 src/Makefile.package.settings（以及 src 目录的 Makefile 体系）中添加PKG_MANYBODY = MANYBODY后执行make yes-manybody。详见 Build package。

7.4 默认值

该 pair style 无默认参数（"none"）。

八、结语与延伸阅读

meam/sw/spline为需要同时体现金属键合与方向性作用的单元素体系（Si、Ti 及其 bcc/hcp 相研究）提供了简洁而精确的势形式：样条化的七函数表示让势函数形状几乎不受解析形式约束，SW 角项又补足了 EAM 家族对角向作用的描述能力。实际使用时只需牢记三条主线——势文件自带截断、mass 单独指定、仅支持单元素——即可快速上手。

进一步阅读建议：

• 势文件仓库：<potentials/Si.b.meam.sw.spline>、<potentials/Ti.meam.sw.spline>
• 源码实现：pair_meam_sw_spline.cpp、pair_meam_sw_spline.h
• 示例目录：examples/PACKAGES/meam_sw_spline（Si 的 bcc/fcc/sc/dc 及对应 relax 版本、Ti 的 α/β/ω/fcc/hexagonal 相脚本）
• 回归测试：atomic-pair-meam_sw_spline.yaml
• 相关命令：pair_coeff、pair_style meam、pair_style meam/spline

参考文献：Lenosky et al.,Modell. Simul. Mater. Sci. Eng.8, 825 (2000)；Stillinger & Weber,Phys. Rev. B31, 5262 (1985)；以及 Nicklas、Fellinger、Park 在俄亥俄州立大学首创的 spline-based MEAM+SW 格式。

• 科研
• 科学计算
• 高性能计算