多核计算机上分子动力学模拟的快速多极子方法并行化

# 多核计算机上分子动力学模拟的快速多极子方法并行化 ## 1. 引言 分子动力学（MD）模拟方法是研究大规模物理/化学系统的传统手段。该方法最初在20世纪50年代提出，但直到70年代中期，随着数字计算机变得强大且价格合理，才开始广泛应用。如今，MD方法仍被广泛用于研究物理/化学系统。 MD本质上是对经典力学（牛顿力学）中的N体问题进行数值求解。对于包含大量粒子N的系统，求解N体问题需要大量时间，并且需要强大的计算机和高效的算法。在解决N体问题时，计算库仑相互作用力是最主要的任务，通常占据MD模拟总计算时间的90 - 95%。 计算力最简单、最准确的算法是直接求和（DS）算法。它使用库仑力公式计算系统中每对粒子之间的相互作用力： \[F_{ij} = k_e\frac{q_iq_j}{||r_{ij}||^2}\cdot\frac{r_{ij}}{||r_{ij}||}\] 其中，\(F_{ij}\) 是位于 \(r_i\) 和 \(r_j\) 处、电荷分别为 \(q_i\) 和 \(q_j\) 的两个粒子之间的库仑力；\(k_e = \frac{1}{4\pi\epsilon_0} = \frac{c_0^2\mu_0}{4\pi} = c_0^210^{-7}H m^{-1}\) 是库仑常数，\(r_{ij} = r_j - r_i\)。在位置 \(r_i\) 处由 \(q_j\) 产生的电势为： \[\Phi_{ij} = k_e\frac{q_j}{||r_{ij}||}\] DS算法的计算复杂度为 \(O(N^2)\)，仅适用于粒子数 \(N\) 最多为 \(10^5\) 的小粒子系统。对于更大的系统，使用DS算法会消耗大量时间。 为了降低N体模拟的力计算成本，人们开发了快速算法，如Barnes - Hut树码（BH）和快速多极子方法（FMM），它们的计算复杂度分别为 \(O(N log N)\) 和 \(O(N)\)。 FMM在许多研究领域有广泛应用，如大规模分子动力学模拟、加速边界元方法、涡方法等。大规模N体或MD模拟，特别是粒子数非常大或不适用周期性边界条件的情况，是FMM应用的典型例子。 由于FMM的广泛适用性，人们进行了许多努力来并行化FMM以实现高性能模拟。目前有几种不同的FMM并行化方法： - 第一种也是最流行的方法是使用消息传递接口（MPI）在分布式内存平台上并行化FMM。 - 第二种是使用专用硬件进行FMM并行化，包括图形处理单元（GPU）、GRAPE（GRAvity piPE）计算机系列等。 - 第三种是使用OpenMP编程模型在多处理器或多核平台上进行并行化。 - 最后一种是结合上述方法。 自多核处理器日益普及以来，已有关于使用OpenMP编程模型并行化FMM的报道。然而，OpenMP的FMM实现数量并不多。主要原因是由于FMM公式复杂，现有的OpenMP实现对于高展开阶数的并行效率中等或较低。此外，随着线程数增加到8或更多，这些实现对于高展开阶数的并行效率会迅速降至40%或更低。 本文介绍了一种在多核计算机上使用OpenMP编程模型实现FMM的方法，以克服现有实现的缺点。开发了一种新的L2L阶段公式和计算过程，简化并加速了FMM的远场力阶段。该方法的主要优点是简单性和并行效率。 ## 2. 快速多极子方法及其变体 ### 2.1 快速多极子方法 FMM是一种 \(O(N)\) 近似算法，用于计算粒子间的力。通过使用多极子和局部展开技术近似力，实现了 \(O(N)\) 的缩放。该算法适用于二维和三维粒子系统。 FMM中力近似的示意图如下： ```mermaid graph LR classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; A(M2M<br>多极子展开):::process --> B(M2L<br>多极子展开到局部展开转换):::process B --> C(L2L<br>局部展开到局部展开转换):::process style A fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; style B fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; style C fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; ``` 一组远处粒子的力通过多极子展开（M2M）进行近似。在观察点，多极子展开转换为局部展开（M2L），然后由观察点周围的每个粒子评估局部展开（L2L）。 使用分层树结构（八叉树）对粒子进行分组。在所有FMM实现中，假设粒子系统位于一个称为根单元的立方体内。根单元被细分为八个相等的子单元，子单元的细分过程继续进行，直到满足某些标准。根单元和子单元形成八叉树，细分过程就是八叉树的构建。当八叉树的层数达到预定义的数量（通常由模拟所需的精度定义）时，八叉树构建停止。 八叉树构建完成后，FMM进入其主要阶段：多极子展开到多极子展开的转换（M2M）、多极子展开到局部展开的转换（M2L）、局部展开到局部展开的转换（L2L），最后是力评估。其中，M2L阶段是最耗时的计算阶段。力评估阶段包括近场力评估和远场力评估两部分。 所有FMM变体都遵循原始FMM的阶段：M2M、M2L、L2L和力评估。唯一的区别是每个变体在这些阶段使用不同的数学公式。以下简要介绍Greengard、Cheng和Rokhlin给出的M2M、M2L和L2L的数学公式： - **球谐函数**： \[Y_m^n(\theta, \varphi) = C\sqrt{\frac{(n - |m|)!}{(n + |m|)!}}P_{|m|}^n(\cos(\theta))e^{im\varphi}\] 其中，\(P_m^n\) 是关联勒让德函数，由Rodrigues公式定义： \[P_m^n(x) = (-1)^m(1 - x^2)^{m/2}\frac{d^m}{dx^m}P_n(x)\] \(P_n(x)\) 表示 \(n\) 阶勒让德多项式。 - **多极子展开**： \[\Phi(X) = \sum_{n = 0}^{\infty}\sum_{m = -n}^{n}\frac{M_m^n}{r^{n + 1}}Y_m^n(\theta, \varphi)\] 其中，\(X_1, X_2, ..., X_n\) 是 \(N\) 个电荷 \(q_1, q_2, ..., q_n\) 的位置，其球坐标分别为 \((\rho_1, \alpha_1, \beta_1), (\rho_2, \alpha_2, \beta_2), ..., (\rho_n, \alpha_n, \beta_n)\)。假设 \(X_1, X_2, ..., X_n\) 位于以原点为中心、半径为 \(a\) 的球体内。点 \(X\) 的球坐标为 \((r, \varphi, \theta)\)。\(M_m^n\) 定义为： \[M_m^n = \sum_{i = 1}^{N}q_i\rho_i^nY_{-m}^n(\alpha_i, \beta_i)\] - **局部展开**： \[\Phi(X) = \sum_{j = 0}^{\infty}\sum_{k = -j}^{j}L_k^jY_k^j(\theta, \varphi)r^j\] 其中： \[L_k^j = \sum_{l = 1}^{N}q_l\frac{Y_{-k}^j(\alpha_l, \beta_l)}{\rho_l^{j + 1}}\] 这里有 \(N\) 个电荷 \(q_1, q_2, ..., q_n\) 位于 \(X_1, X_2, ..., X_n\) 处，其球坐标分别为 \((\rho_1, \alpha_1, \beta_1), (\rho_2, \alpha_2, \beta_2), ..., (\rho_n, \alp