LAMMPS 并行算法深度解析：域分解、Ghost 原子通信、邻居表与并行 FFT 的实现原理

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

LAMMPS 采用基于 MPI 的分布式内存并行架构，通过空间域分解（domain decomposition）+ ghost 原子通信 + Verlet 邻居表三级机制，实现对大规模原子体系的强扩展与弱扩展并行。本文以 LAMMPS 开发者手册中的并行算法章节为核心，结合仓库中src/下的实际源码（comm_brick、comm_tiled、邻居表、KSPACE/FFT 与 OPENMP 包），系统讲解域分解方案（brick/tiled）、forward/reverse 通信模式、邻居表与 bin/stencil 构建算法、PPPM 长程静电力的并行 FFT，以及 MPI 内层的 OpenMP 线程并行设计，帮助开发者完整理解 LAMMPS 并行的底层原理并指导负载优化实践。

一、整体并行架构：强扩展与弱扩展

LAMMPS 的设计目标是支持 MPI 标准下的分布式内存并行，通过域分解把模拟盒拆分为若干互不重叠的子域（subdomain），每个 MPI 进程"拥有"（own）一个子域及其中的原子。其并行化追求两类扩展性：

• 强扩展（strong scaling）：同一体系下，增加进程数能近似按比例加速；
• 弱扩展（weak scaling）：体系规模扩大时，计算成本与体系规模成正比。

在此之上，还可以在同一 MPI 子域内部使用 GPU 或 OpenMP 做第二层并行。手册章节为清晰起见多用 2d 示意图，但所述算法在 2d 和 3d 情形完全通用。

注意：该章节文字与配图改编自 LAMMPS 正式发表论文中"并行算法"一节，见文档中的 论文引用。

二、域分解（Partitioning）：brick 与 tiled 两种模式

2.1 模拟盒与两种分解模式

LAMMPS 的模拟盒是一个 3d（或 2d）体积，可以是**正交（orthogonal）或单斜/三斜（triclinic）**形状：

• 正交盒：盒边沿x、y、z笛卡尔轴对齐，各面均为矩形；
• 三斜盒：盒边沿三个任意向量对齐，各面为平行四边形。

每个维度上的盒面可以是周期边界的，也可以是非周期的（fixed 固定边界：原子移出即删除；shrink-wrapped：盒面随原子包络连续收缩调整）。

分解方式由comm_style命令选定，只有两种：

正交盒与三斜盒都可以在模拟中连续变形（例如压缩固体、剪切液体），处理器子域会随之变形。源码层面，src/comm.h 中的Comm基类明确定义了BRICK（6 向 stencil 通信）与TILED（不规则瓦片通信）两种 style，以及LAYOUT_UNIFORM / LAYOUT_NONUNIFORM / LAYOUT_TILED三种布局；src/comm_brick.cpp 与 src/comm_tiled.cpp 分别实现两种分解。处理器网格的划分逻辑则位于 src/domain.cpp（Domain类负责依据体系密度比与进程数优化procgrid[3]）。

2.2 非均匀密度下的负载不均衡与三种缓解策略

密度不均匀的模型中，默认 brick 分解会导致各进程粒子数不均衡，整体模拟速率受最慢（负载最大）进程限制。手册给出了三种官方缓解手段：

1. 调整处理器网格：processor grid 默认基于模拟盒体积、并尝试优化子域的"体积/表面积"比，可用processors命令手动指定Nx Ny Nz；
2. 平移分解平面：balance命令可以平移定义子域大小的并行平面，常与更优的处理器网格组合使用；
3. 递归二分（RCB）+ tiled 通信：balance的 RCB 模式结合comm_style tiled，可产生每个子域粒子数相等的划分。

这些实现在源码中一一对应：src/balance.cpp 及其imbalance_*系列文件实现了各种平衡算法（包括 RCB），src/comm_tiled.cpp 中包含递归遍历 RCB 切割树（rcbinfo数组记录每个进程的切割维度与切割比例cutfrac）来判定原子归属、计算子域间重叠的逻辑。

手册用一个含真空的 2d 体系（12 个 MPI rank）演示了策略差异：默认分解下多个 rank 分不到原子（严重失衡）；强制1x12x1网格后每个 rank 都有计算量但子域内原子数仍不均，且薄片形状增大通信量；2x6x1网格并平移分解平面后负载与通信进一步改善；RCB 递归二分给出最均衡的划分（原子颜色从绿到红表示子域负载优劣）。

三、Ghost 原子通信：forward 与 reverse 两种模式

3.1 为什么需要 ghost 原子

按选定的分解方案，所有 per-atom 数据都分布在 MPI 进程之间——这正是 LAMMPS 能用足够多的 MPI 进程处理超大体系的前提。但计算短程作用力时，每个进程不仅需要访问自己拥有的原子，还需要邻域子域中的原子信息，即 LAMMPS 所说的ghost 原子：邻居进程原子的副本，保存最多到通信截断距离（communication cutoff）内的 per-atom 数据。

这一机制同时也是周期性边界条件的实现方式：位于周期边界另一侧、距离在截断之内的原子同样存为 ghost 原子，取自子域的周期复制（串行运行时即取自本进程自己）。由此带来的一个重要性质是：LAMMPS 的力计算不受最小镜像约定（minimum image convention）限制，截断距离可以大于模拟域尺寸的一半。

3.2 两阶段（3d 为三阶段）交换模式

由于 ghost 数据必须在每个 MD 时间步或最小化步更新，通信模式必须高效。2d 情形下通信分两阶段（3d 分三阶段）进行，规则与不规则分解均适用：

• 规则（brick）情形：先沿x方向交换，再沿y方向，处理器子域网格中每个进程有 4 个邻居。x阶段中，rank 1 和 2 把红色阴影区域中的 owned 原子发给 rank 0（反之亦然）；y阶段中，rank 3 和 4 把蓝色阴影区域的原子发给 rank 0——这些原子包含它们在x阶段收到的 ghost 原子。rank 0 由此集齐全部 ghost 原子；图中实心蓝色角区里的原子会被通信两次后才到达 rank 0。
• 不规则（tiled）情形：两阶段逻辑相似，但一个处理器在每个方向上可能有多个邻居。x阶段中 rank 1、2、3 把红色区域原子发给 rank 0，且只包含位于 rank 0 子域y上下边界之间的原子；y阶段中 rank 4、5、6 把蓝色区域原子发给 rank 0，其中可能包含x阶段收到的 ghost 原子，但仅限 rank 0 填充其 ±y方向扩展 ghost 区域所必需的。因此 rank 5、6 之间互相收到的 ghost 原子不会转发给 rank 0。要点是：尽管不规则分解的通信模式更复杂，仍然只需两阶段（3d 三阶段）、仅与相邻处理器交换即可完成。

3.3 forward / reverse 通信的定义与实现细节

• forward 通信：owned 原子的属性发给邻居进程、成为其 ghost 原子的属性。在每个原子迁移、邻居表重建的时间步，每个处理器会列出"自己是哪些邻居进程 ghost 原子"的 owned 原子清单；此后直到下次 reneighboring，都用这些清单把 per-atom 坐标等数据打包进消息缓冲。
• reverse 通信：计算出的 ghost 原子属性送回原子所属进程，典型用途是把 ghost 原子上的部分力累加回 owned 原子的完整力。两阶段顺序与 forward 相反，使用同一套原子清单打包/解包力数据；缓冲解包时 ghost 力直接加到 owned 原子力上。图中四个蓝色角区原子的力在两个阶段各发送/接收/累加一次，owned 处理器最终才拿到完整力。

这两个操作在 LAMMPS 内部被广泛复用：坐标与力交换之外，manybody 势用它共享 per-atom 中间量、刚体积分器用它让 body 内各原子访问 body 属性。

手册还列出若干底层细节，均能在源码中印证：

• 与不同处理器交换时使用MPI_Send()/MPI_Irecv()；若处理器"与自己交换"（如单进程跨周期边界造 ghost），只做 pack/unpack。src/comm_brick.cpp 的forward_comm()与reverse_comm()正是这种"先MPI_Irecv发异步接收、再MPI_Send发数据"的经典实现。
• 跨周期边界发送 owned 原子坐标时，收发过程中直接加减周期盒长，因此建邻居表与算力时无需再做最小镜像。
• ghost 交换截断通常等于邻居截断，但可设得更大：如多原子刚体的半直径、扭结作用力的 3 倍以上拉伸键长，甚至可以超过周期盒尺寸。规则模式下，若截断超出某方向邻居子域范围，会沿同一方向做多次交换（第一次只发 owned 原子，后续各次发送之前收到的 ghost）；不规则模式则检测扩展 ghost 子域与其他处理器在每个方向上的全部重叠。
• Comm基类中的cutghost[3]、cutghostuser以及SINGLE/MULTI多集合截断模式（src/comm.h第 31–42 行）对应了"单截断 / 多集合截断"两种通信模式；forward_comm()/reverse_comm()为纯虚函数并被Pair、Bond、Fix、Compute、Dump等重载（src/comm.h），供各组件复用同一通信基础设施。

四、邻居表（Neighbor Lists）：Verlet 列表、bin 与 stencil

4.1 Verlet 列表与 skin 距离

为高效算力，每个处理器构建 Verlet 风格邻居表，枚举所有i,j对（i= owned，j= owned 或 ghost）中距离小于邻居截断的原子对。实现要点：

• 邻居表存储在多页（multiple-page）数据结构中：每页是一段连续内存，为多个i原子存放邻居j的向量，可整块分配/释放。源码对应 src/my_page.cpp 中的MyPage页分配器与 src/neigh_list.cpp 中的NeighList类。邻居表通常是 LAMMPS 中内存占用最大的数据结构。
• 邻居表每隔若干步（从零）重建一次，其间每步重复使用。邻居截断距离为 $R_n = R_f + \Delta_s$：$R_f$ 是势函数定义的最大力截断，$\Delta_s$ 是skin 距离，允许假设原子相邻步间不会移动太远而多步复用列表。
• 通常当任一原子自上次重建以来移动了半 skin 距离时触发重建；相关选项由neigh_modify命令调整，原子排序则由atom_modify控制。

4.2 重建步：原子迁移、周期性包裹与空间排序

执行 reneighboring 的时间步上：

1. 移出所属子域的原子先经通信迁移到新进程；周期性边界条件也仅在这些步执行，保证每个原子被分配到正确进程；
2. 迁移后，owned 原子存入连续向量；
3. 每隔若干步对 owned 原子做空间排序（spatial sort）：先用 bin 分箱，再重排使同 bin 原子在向量中相邻，以改善算力与建表时的缓存效率。

4.3 bin 网格与 stencil 搜索

为线性时间构建本地邻居表，模拟域上叠加规则 3d（或 2d）bin 网格。每个处理器存储覆盖"子域 + 邻居截断 $R_n$"的 bin 集合。bin无需对齐处理器边界：每个维度的整数 bin 数是按整个模拟盒尺寸整除适配的。

最常见的构建是半列表（half list）：每对i,j只存一次。为此使用围绕中心 bin 的预计算stencil——一组x,y,z整数偏移，列出足够近、可能包含任何距中心 bin 内原子i在 $R_n$ 范围内的邻居j的 bin。半列表的 stencil 可以是不对称的。

手册给出的 stencil 数值（bin 尺寸为 $\frac{1}{2}R_n$ 时）：

三斜 stencil 多出的 bin 是因为 bin 平铺的是三斜域的包围盒，本身相对于模拟盒并不周期，故判定哪些i,j对入表的逻辑略有调整。

构建算法：先遍历 owned + ghost 原子并分配 bin（用整型向量在每个 bin 内建原子索引链表），再三重嵌套循环——owned 原子i→i所在 bin 的 stencil → 每个 stencil bin 内的j原子（含 ghost）——若 $r_{ij} < R_n$ 则把j加入i的邻居向量。

4.4 邻居表构建的扩展选项

手册列出的构建选项与源码/命令对应关系：

• bin 尺寸可调：$\frac{1}{2}R_n$ 在多数典型情形最优——bin 更小则循环开销增大，bin 更大则距离计算变多。bin 更小时 2d stencil 更接近半圆（3d 半球），方形角上的 bin 因距离过远被剔除；
• 多列表与多种 stencil：不同截断（算力 vs 如 RDF 诊断）、r-RESPA 按距离分层的子集、"full" 列表（i,j存两次，用更大的对称 stencil）、部分枚举 ghost 邻居的列表——后两者被各类 manybody 势使用；入表判据除距离外还可依赖两粒子直径（如多分散颗粒）；
• pair_style hybrid：需为主邻居表按子势生成子列表，每个子列表只含相应原子子集间的i,j对；
• 异构截断体系：如溶胶-胶体体系（colloid/colloid、colloid/solvent、solvent/solvent 截断差异巨大）、粒径比可达 10–100 倍的多分散颗粒。LAMMPS 提供高效算法：对不同截断使用不同 stencil 并裁剪出仅跨越截断球面的 bin；还有更新的多 stencil + 多 bin 尺寸方法，可高效处理任意粒径比体系；
• neighbor nsq：小/稀疏体系或兜底方法，无分箱、子域内全 $O(N^2)$ 双重循环；
• newton pair开关的影响：newton pair off时半列表包含所有j为 ghost 的对；newton pair on时仅当j的z（或z相等时y、再相等时x）坐标更大才入表。对均匀致密体系，这使 owned 原子固定从同一方向扇区取邻居，邻居表长度更一致、降低负载不均衡风险。

五、长程相互作用的并行实现：PPPM 与 3d 并行 FFT

5.1 PPPM 的五阶段流程

带电体系的长程库仑作用由 FFT 基 PPPM 方法（kspace_style pppm及其变体）计算：库仑作用拆为短程部分（实空间、邻居表内成对循环）与长程部分（倒易空间、kspace 风格计算）。PPPM 在模拟盒上叠加 3d 规则 FFT 网格，流程为：

1. 每个原子的点电荷插值到邻近网格点；
2. 前向 3d FFT；
3. 倒易空间卷积；
4. 一个或多个逆 3d FFT；
5. 由原子附近网格点的电场值插值回原子，得到力。

5.2 网格的域分解与 grid ghost

无论哪种空间分解，每个处理器都拥有其子域内的砖形（brick）网格点集合。两次插值都用到围绕每个原子的网格点 stencil，为此每个处理器额外保存包围其砖块的若干层ghost 网格点。网格点值的 forward/reverse 通信与原子数据通信方式类似：owned 网格点上的电场发给邻居成为其 ghost 点值；ghost 点上的电荷值发送并累加回 owned 点。

三斜盒情形：FFT 网格平面平行于盒面，但电荷/电场的网格映射在约化坐标中进行（把倾斜盒概念上视为单位立方体），stencil 与 FFT 操作不变；不过达到给定精度所需的 FFT 网格比正交盒更大。

5.3 -pencil 分解的并行 FFT

并行 3d FFT 的通信量相对其计算量很大。3d FFT 实现为沿x、y、z三个方向的 1d FFT 序列，因此网格不能像原子那样按 3 维分解，只能按 1 维（平面）或 2 维（铅笔）分解，步骤之间需要转置，使每个 MPI 进程"拥有"的网格部分恰好完整覆盖其待做 1d FFT 的方向。LAMMPS 采用图中所示的pencil 分解：

1. 初始（最左）：每处理器拥有子域内的同色网格砖；
2. brick-to-pencil 通信：转为x方向 1d 铅笔布局（左二），随后用配置的 FFT 库对每条完整拥有的铅笔做 1d FFT；
3. pencil-to-pencil 通信：转为y铅笔（相当于转置x、y维），做y方向 1d FFT；
4. 最后一次y→z铅笔转置 +z方向 1d FFT，完成前向 FFT；数据保持z-pencil 布局供卷积；
5. 逆 FFT 按相反顺序执行 1d FFT 与通信，最终网格值布局与初始砖布局相同。

每次 FFT 内部通信（brick-to-pencil、pencil-to-pencil、pencil-to-brick）都是把一种tiling（把 3d 网格的小砖形子集分配给各处理器的覆盖方案）转换为另一种。LAMMPS 配套的fftMPI库支持任意 tiling 定义，因此不规则域分解可直接使用；tiling 转换用点对点通信实现——每个处理器只向少数处理器发送数据（初始 tiling 中每个 tile 只与目标 tiling 的少数几个 tile 重叠）。

虽然理论上可用单次MPI_Alltoall()做全体 $P$ 个进程集合通信，但通常慢得多。对图中 brick/pencil tiling，甚至不需要全体集合通信：$8^3$ 网格（512 网格单元）分到 64 处理器、每处理器拥有一块 2×2×2 砖时，brick-to-pencil 只需在4 个处理器的小组内做集合通信（图中 4 种颜色）；一般地，可把 $P$ 个处理器分为 $P^{2/3}$ 个小组、每组 $P^{1/3}$ 个处理器，各小组内部做集合通信。pencil-to-pencil 则分为 $P^{1/2}$ 个小组、每组 $P^{1/2}$ 个处理器，如图示y⇒z情形中 8 处理器小组拥有前部yz平面并在组内转置。

默认使用点对点通信；kspace_modify collective yes可切换为分区集合通信模式，此时代码检测不相交子组大小，把单个 $P$ 大小的 communicator 划分为多个小 communicator 分别调用集合通信。文档引用了 Argonne 国家实验室 IBM Blue Gene/Q 上的测试：大进程数下分区集合通信的 FFT 性能显著优于纯点对点或全体 $P$ 进程集合通信。

相关 FFT 实现位于 KSPACE 包，如 src/KSPACE/fft3d.cpp、src/KSPACE/fft3d_wrap.cpp，KOKKOS 包另有 src/KOKKOS/fft3d_kokkos.cpp。

手册补充的 FFT 细节：

• fftMPI允许每个网格维度是小素因子（2、3、5）的倍数，支持任意处理器数；brick/pencil 分解不一定对齐，但两种通信的通信子组规模仍按 $O(P^{1/3})$ 与 $O(P^{1/2})$ 缩放；
• 为 1d FFT 效率，转置操作同时重排 3d 数据使不同方向在内存中连续，重排可并入 MPI 缓冲的 pack/unpack；
• 大体系、多进程下并行 FFT 的主导成本常是通信而非 1d FFT 计算（后者按每方向网格点数 $N$ 计为 $N\log N$），因为网格只能 2d 铅笔分解，而原子数据（短程力/能计算）可高效 3d 分解。两种对策：
  • MPI + OpenMP 混合并行（package omp）：减少参与 FFT 通信的 MPI 进程数，MPI 域内用 OpenMP 保持算力。部分计算并行效率略降可能小于 FFT 通信开销；
  • 多分区 +run_style verlet/split：启动多分区（multi-partition）计算，用整数 "1:p" 映射把p个原子子域映射到 1 个 FFT 网格子域，非键合/键合力与能在较大分区计算，PPPM kspace 计算并发地在较小专用分区执行。

此外：LAMMPS 还有基于 PPPM 的偶极、色散（Lennard-Jones）长程作用求解器，可与点电荷库仑长程作用联用；GridComm类在模拟域上叠加规则网格并按与处理器子域一致的方式分区，提供 owned/ghost 网格点值的 forward/reverse 通信方法，被 PPPM（FFT 网格）与 MSM（从细到粗的网格级联库仑力算法）使用，也适用于连续场与粒子耦合的模型，例如两温模型（TTM）中fix ttm/grid用GridComm在分布式网格上做有限差分热方程（替代使用复制网格的fix ttm）。

六、第二层并行：OPENMP 包的设计与数据竞争规避

INTEL、KOKKOS、OPENMP 三个包均提供 OpenMP 线程并行，主要把循环分布在本地数据上，与 MPI 空间域分解构成正交的两层并行策略。手册以结构最简单的OPENMP 包为例讲解（INTEL/KOKKOS 更复杂，支持向量化、加速器与双精度/单精度/混合精度选择；其写法示例见开发者手册中"为 OPENMP 包编写 style"一节）。

6.1 引入 OpenMP 的理由

• 资源利用效率：多核机器上 MPI 任务常受通信带宽限制（节点间通信共享同一网络链路）；减少 MPI 任务、增加 OpenMP 线程可降低通信开销与内存争用。代价是节点内内存带宽需求更高——域分解天然利于 CPU 缓存局部性，而多线程实现为规避数据竞争使用 per-thread 复制存储，需要额外的规约（reduction），线程数少时效率最高；
• PPPM 扩展性：长程静电求解器在高节点数下受 3d FFT all-to-all 通信限制（子域变小后计算/通信比恶化），用 OpenMP 后求解器可在 MPI 任务子集上运行，改善整体扩展；
• 并行粒度：OpenMP 提供节点内按粒子/邻居的更细粒度并行；
• 负载均衡：域分解隐式假设体系均匀；非均匀体系中域分解的负载均衡手段也有局限。改用 OpenMP 后每个 MPI rank 的子域更大，负载均衡通常更有效；
• 高性能计算实践：混合 MPI+OpenMP 常是大型 HPC 集群/超算获得最优性能的关键。

6.2 包的设计：继承 + 最小修改 + 辅助类

OPENMP 包的关键设计决策是保持源码改动最小，便于维护并与非线程版标准实现保持同步：

• 继承：OPENMP 版本从标准版本派生（如PairLJCutOMP继承自PairLJCut），只覆写多线程化或需修改的方法（沿用 OPT 包的做法）。仓库中 src/OPENMP/ 目录包含上百个此类派生类，如 src/OPENMP/pair_lj_cut_omp.h 等；
• 最小修改：修改后的代码结构与标准版尽量保持一致，便于逐行对照；
• 辅助类：多线程支持功能拆分到三个类——
  • ThrOMP：提供"_thr" 变体的 tally 等多线程感知功能（多重继承），见 src/OPENMP/thr_omp.h；
  • ThrData：管理 per-thread 数据结构，避免 false sharing 拖慢，见 src/OPENMP/thr_data.h；
  • FixOMP：管理由package omp命令激活的全局多线程状态与设置，见 src/OPENMP/fix_omp.h。

6.3 规避数据竞争的三条路线与 LAMMPS 的选择

MD 代码线程并行的核心问题是更新力、能、应力、力矩等 per-atom 累积量时的数据竞争：交互计算总是涉及多个原子，多线程会同时更新同一原子的数据。不做大幅重写的前提下有三条路线：

1. 锁或原子操作：保证同一时刻只有一个线程更新，适合冲突很少的场景；
2. 每个"拥有"元组中原子的线程都计算全部 per-atom 属性，只存 owned 原子的数据，线程多时最有效，GPU 加速中常用；
3. per-thread 数据副本：各线程独立更新自己的副本，最后用 reduction 合并回常规 per-atom 存储，需要额外的初始化/清零/合并步骤。

OPENMP 包选择复制的 per-thread 数据结构，因为：单线程时性能无损（区别于另两方案）；代码可维护且与非线程版相似；且对 2–8 线程这种 OpenMP 的典型场景效率最高。

6.4 循环调度与邻居表并行

• 循环调度：多线程并行通过把（外层）循环静态分配到各线程实现，通常是处理邻居表i,j对或 bond/angle/dihedral/improper 列表时对本地原子i的循环。由于均匀体系中每原子邻居数相近，线程间显著负载不均衡不常见。
• 邻居表构建并行化：除力计算外，邻居表生成本身也被并行化——每个线程处理一块不同的 owned 原子，并用MyPage页分配器实例管理自己的一组邻居"页"。

七、实践要点小结

从源码结构看，上述机制的落点分别是：分解与通信在 src/comm_brick.cpp 和 src/comm_tiled.cpp，进程网格在 src/domain.cpp，负载均衡在 src/balance.cpp，邻居表与页分配在 src/neigh_list.cpp 与 src/my_page.cpp，并行 FFT 在 src/KSPACE/fft3d.cpp 系列文件，OpenMP 层在 src/OPENMP/ 目录。结合这些文件阅读本文各节，可以完整掌握 LAMMPS 并行算法从域分解、原子通信、邻居表到长程作用与多线程实现的全链路原理。

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