LAMMPS compute composition/atom：逐原子局部组成向量的用法、实现原理与测试验证
LAMMPS compute composition/atom：逐原子局部组成向量的用法、实现原理与测试验证
本篇指南基于 LAMMPS 官方文档 doc/src/compute_composition_atom.rst 及其配套源码，系统讲解 compute composition/atom 命令：它如何为每个原子在给定截断半径内统计邻域中各原子类型的占比（局部组成向量），cutoff 参数与 pair style、comm_modify 之间的联动规则，输出数组的列含义，以及 OPENMP / KOKKOS 加速变体的实现差异。读完后你可以直接在合金、共混、溶质-溶剂等模拟场景中部署逐原子局部成分分析，并能理解其邻域统计的实现边界与常见陷阱。
命令语法与参数
compute composition/atom 属于 per-atom 型 compute，其命令格式为：
compute ID group-ID composition/atom keyword values ...
其中 ID 与 group-ID 的含义见 compute 命令文档。该命令只接受一个可选关键字：
| 关键字 | 取值 | 含义 |
|---|---|---|
cutoff |
数值（> 0） | 邻域搜索所用的距离截断半径 |
官方文档给出的最小示例：
compute 1 all composition/atom
compute 1 all composition/atom cutoff 9.0
comm_modify cutoff 9.0
第二个示例特意搭配了 comm_modify cutoff 9.0，原因在于：当指定的 cutoff 大于 pair style 截断半径加上 neighbor skin（或根本没有定义 pair style）时，必须用 comm_modify cutoff 把通信截断同步设置为相同值，否则跨处理器的 ghost 原子不足以覆盖搜索半径。这一约束并非只是文档约定，而是在源码初始化阶段强制校验的——ComputeCompositionAtom::init() 中会计算 cutghost = MAX(pair->cutforce + skin, comm->cutghostuser)，若 cutoff > cutghost 则直接报错 Compute composition/atom cutoff exceeds ghost atom range - use comm_modify cutoff command（源码第 91-94 行）。
参数解析逻辑同样值得注意（构造函数第 39-63 行）：
- 命令总参数个数必须在 3～5 之间（
ID group-ID [cutoff value]），否则报Illegal compute composition/atom command； cutoff值必须为正数，且必须能解析为合法数值；- 未指定
cutoff时，内部将cutoff置为 0，留待init()中取 pair style 的cutforce作为默认值； - 构造函数中即按
size_peratom_cols = 1 + ntypes预定好输出列数（第 59-60 行），ntypes即atom_style定义的原子类型总数。
局部组成的物理定义
该 compute 于 LAMMPS 21Nov2023 版本加入。对每个中心原子，设其在截断球内有 M 个邻居，则“局部组成”定义为：某一原子类型的原子个数（含中心原子本身）除以 (M+1)。对给定中心原子，所有类型占比之和恒等于 1。
文档同时给出一个重要的使用提示：本 compute 使用的是 pair_coeff 中分配的原子类型（atom type）数，而不是化学元素（chemical species）。如果一个势函数把两种元素（例如 Cu 和 Ni）分配给了四个不同的原子类型（如 'Cu Cu Ni Ni'），计算结果就是 4 个分量的占比向量。若你希望得到按元素合并的结果，需要额外一步聚合：可用 compute reduce、原子型 variable 在 LAMMPS 内部完成，也可以作为后处理步骤离线合并。
这一“按类型累加、最后归一化”的定义可以直接在源码中对应到实现上：compute_peratom() 对 group 内每个原子先令 count = 1（中心原子自身）并把 result<a href="https://link.gitcode.com/i/7fe332188c842bd53fd65c477935fe23" target="_blank">i][itype] 加 1，随后遍历邻域列表中满足 rsq < cutsq 的每个邻居 j，对 result[i][jtype] 累加并 count++；最后写 result[i][0] = count，并以 lfac = 1.0 / count 将第 2 列起的各类型计数统一乘为分数（[第 178-206 行）。
输出信息：1+N 列的 per-atom 数组
该 compute 产生一个每原子 1 + N 列的数组，N 为原子类型数：
- 第 1 列：参与组成统计的原子总数（含中心原子）；
- 第 2 … N+1 列：每种原子类型在截断球内的占比（0～1 的分数，各行之和为 1）。
这些值可以被任何接受 compute 原子值的命令消费，例如 dump custom、fix ave/time 等，LAMMPS 输出选项的总览见 Howto output。不属于 group-ID 的原子输出整行零值——源码中先对整个 result 数组 memset 清零，再只对 mask<a href="https://link.gitcode.com/i/76ccc4589fe8f73255b37f2b473c7c9d" target="_blank">i] & groupbit 为真的原子写数（[第 151-168 行），这一细节在官方单元测试数据中也得到了印证（见下文“单元测试”一节，非 group 原子的 6 列输出全为 0）。
cutoff 的默认值与邻域列表复用策略
关于 cutoff 的默认行为，文档规定为：默认取 pair style 指定的截断半径；若没有定义任何 pair style，则必须通过该关键字指定 cutoff。这一点由 init() 的前几行 强制保证：if (!force->pair && cutoff == 0.0) 直接报错 Compute composition/atom requires a cutoff be specified or a pair style be defined。
从源码结构看，该 compute 还包含一层“邻域列表复用”优化（init() 第 97-112 行）：
- 若
cutoff == 0，则以pair->cutforce作为实际截断； - 若最终
cutoff <= pair->cutforce + skin，令cutflag = 0，即可以复用 pair style 已经建好的完整邻域列表；否则cutflag = 1； - 向邻居系统提交请求
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL)——这是一个“完整列表 + 按需（occasional）构建”的请求；仅当cutflag = 1时才通过req->set_cutoff_fixed(cutoff)用自定义截断单独建表。
同时，init() 明确拒绝 neighbor style 为 multi 的构建：Compute composition/atom requires neighbor style 'bin' or 'nsq'（第 105-106 行），与文档“Restrictions”一节的要求一致，即该 compute 要求 neighbor 使用 bin 或 nsq。
另一个来自文档的重要性能提示是：计算本量所需的邻域列表会在每次执行计算时重建（即每次对该量做快照，通常是 dump 触发时）。因此过于频繁地计算/转储该量是低效的。这与源码中每次 compute_peratom() 都调用 neighbor->build_one(list)（第 144 行）相吻合——build_one 会按需复制或重新构建该请求对应的列表。实践中建议只在 dump custom 的间隔内按需采样该组成量，而不是每个时间步都计算。
bonded 体系中的 special_bonds 陷阱
文档“Description”中的第二个 note 给出了一个容易被忽视的边界情况：若体系存在键合相互作用，special_bonds 命令的设置会使同一 bond/angle/dihedral 内的原子对从相互作用中剔除。由于 special_bonds 的默认设置即如此，这些原子对不会出现在邻域列表中；而 composition/atom 依赖邻域列表工作，因此这些被剔除的原子对也不会被计入局部组成。
文档给出的规避方案是：写出 dump 文件，再用 rerun 命令对 dump 中的快照重新计算组成，并在 rerun 输入脚本中把 special_bonds 设为包含邻域列表中的所有键合对。官方单元测试正是采用了类似思路：其 in.fourmol 输入脚本使用 pair_style zero 8.0 等 zero 风格（不产生真实相互作用），并显式设置 special_bonds lj/coul ${bond_factor} ${angle_factor} ${dihedral_factor}，保证邻域列表的构成完全由几何截断决定，便于对组成向量做确定性校验。
加速变体：composition/atom/omp 与 composition/atom/kk
文档顶部声明了加速变体 composition/atom/kk 和 composition/atom/omp。按照 LAMMPS 加速器包的通用规则（见文档尾部包含的 accel_styles.rst）：带 gpu、intel、kk、omp、opt 后缀的风格在功能上与无后缀版本相同，只是针对可用硬件做了加速；带后缀的风格接收相同参数、应产生相同结果（除舍入与精度问题外）。它们分别属于 GPU、INTEL、KOKKOS、OPENMP、OPT 包，仅在 LAMMPS 以相应包构建时才可用（构建方式见 Build package）。你既可以在脚本中显式写后缀，也可以在调用 LAMMPS 时通过 -suffix 命令行开关（见 Run options）或在输入脚本中用 suffix 命令自动切换。
两种加速实现的具体做法：
- OPENMP 变体（src/OPENMP/compute_composition_atom_omp.cpp）：仅把逐原子的外层循环用
#pragma omp parallel for schedule(static)并行化（第 76-79 行）。源码头部的注释明确说明：每个原子只写自己的输出行result[i]，因此无论线程数或调度方式如何，结果都与串行版本逐位一致（bit-identical）。 - KOKKOS 变体（src/KOKKOS/compute_composition_atom_kokkos.h 注册了
composition/atom/kk、composition/atom/kk/device、composition/atom/kk/host三个风格名）：在 init() 中根据设备类型调整邻居列表请求的kokkos_host/kokkos_device标志；计算主体改用Kokkos::parallel_for("ComputeLocalComp", ...)将同一套“计数—归一化”逻辑映射到设备上执行（operator() 第 115-165 行）。注意设备端累加使用KK_FLOAT（单精度）类型，这与文档“除舍入与精度问题外结果相同”的表述相对应：GPU 单精度构建下输出可能存在小的数值偏差。
单元测试中的完整算例
仓库自带的单元测试为验证输出格式提供了最直接的参照：unittest/force-styles/tests/compute-composition_atom.yaml。该测试的前置条件是 atom full + compute composition/atom，在 in.fourmol 脚本末尾追加：
compute test solute composition/atom cutoff 4.0
其中 in.fourmol 使用 pair_style zero 8.0（配合 coeffs.fourmol 中定义的 lj/cut 8.0 势与 5 种原子类型），读取 29 个原子的 data.fourmol。由于指定了 cutoff 4.0（小于 pair 截断 8.0 + skin），无需额外 comm_modify。参考数据展示了“1+N = 6 列”输出的具体形态，例如：
1.0000000000000000e+00 1.1000000000000000e+01 1.8181818181818182e-01 5.4545454545454541e-01 1.8181818181818182e-01 9.0909090909090912e-02 0.0000000000000000e+00
即：原子 1 在 4.0 截断球内共 11 个原子（第 1 列），类型 2 占比约 0.545（=6/11）、类型 1 和类型 3 各 2/11、类型 4 为 1/11，各列占比之和为 1。数据中后 12 个原子（18～29 号，属于 solvent 而非 solute group）输出全零行，直观演示了“非 group 原子只报告零值”的行为。
构建限制与相关命令
- 包依赖：该 compute 属于 EXTRA-COMPUTE 包，只有以该包构建的 LAMMPS 才可用（见 Build package）。
- 邻域风格：要求
neighbor为bin或nsq（源码中显式拒绝multi，见 init() 第 105-106 行）。 - 默认值：
cutoff默认为 pair style 截断半径（pair->cutforce）。 - 相关命令：comm_modify（大截断场景必须同步）、special_bonds（键合体系的邻域剔除问题）、rerun（快照重算方案）、pair_coeff（决定类型计数口径）。
小结
compute composition/atom 是一个典型的“基于邻域列表的逐原子统计”实现：文档层面定义了截断球内“类型数/(M+1)”的组成口径与 1+N 列输出；源码层面（compute_composition_atom.h / compute_composition_atom.cpp）则补充了文档未展开的工程细节——cutoff 与 ghost 范围的强校验、pair 列表复用与 REQ_FULL | REQ_OCCASIONAL 按需建表、非 group 原子输出零行、以及每次计算都会触发 build_one 重建列表的性能特征。使用它时把握三条要点即可：自定义截断超过 pair 截断 + skin 时同步 comm_modify cutoff；在键合体系中意识到 special_bonds 对邻域的剔除效应；以及避免对该量做过于高频的 dump/计算。