C++ 分子动力学模拟代码

Question

我正在开展一个模拟气体硬球模型的项目。 （类似于理想气体模型。）

我已经编写了整个项目，并且它正在运行。为了让您了解我所做的事情，有一个循环执行以下操作：（伪代码）

Get_Next_Collision(); // Figure out when the next collision will occur
Step_Time_Forwards(); // Step to time of collision
Process_Collision(); // Process collision between 2 particles
(Repeat)

对于大量粒子（比如 N 个粒子），必须进行 O(N*N) 检查以确定下一次碰撞发生的时间。遵循上述程序显然是低效的，因为在绝大多数情况下，粒子对之间的碰撞不受其他地方的碰撞处理的影响。因此，希望有某种形式的优先级队列来存储每个粒子的下一个事件。 （实际上，由于一次碰撞涉及 2 个粒子，因此只会存储一半的事件，因为如果 A 与 B 发生碰撞，那么 B 也会与 A 发生碰撞，而且是同时发生的。）

我发现编写这样一个事件/冲突优先级队列很困难。

我想知道是否有任何已编写的分子动力学模拟器，我可以去查看源代码以了解如何实现这样的优先级队列。

通过谷歌搜索，我很清楚已经编写了许多 MD 程序，但其中许多程序要么过于复杂，要么不适合。

这可能是因为它们具有巨大的功能，包括生成可视化的能力或计算粒子之间的相互作用力的模拟的能力等。

有些模拟器不适合，因为它们为不同的模型进行计算，即：不是能量守恒的硬球模型，具有弹性碰撞。例如，与势能相互作用的粒子或非球形粒子。

我曾尝试查看 LAMMPS 的源代码，但它非常庞大，我很难理解它。

我希望这是关于我正在尝试做的事情的足够信息。如果没有，我可能会添加更多信息。

Answer 1

位置感知系统的基本版本可能如下所示：

1. 将宇宙划分为一个立方体网格（其中每个立方体都有边 A，体积为 A^3），其中每个立方体都足够大，但比系统的总体积要小得多。每个网格立方体进一步分为 4 个子立方体，其粒子理论上可以提供给相邻的立方体（并借给计算）。
2. 每个网格立方体都会记录其中包含的粒子，并知道其相邻网格立方体包含的粒子。
3. 将粒子的可观测宇宙定义为半径为（网格维度/2）。定义timestep=(griddim/2) / max_speed。这假设来自最多四个相邻网格立方体的粒子理论上可以在该时间段内相互作用。
4. 对于每个网格立方体中的每个粒子，运行传统的碰撞检测算法（使用 mini_timestep < timestep，检查每个粒子是否与其可观察宇宙中的其他粒子可能发生碰撞。将碰撞存储到任何按时间排序的结构，甚至只是一个数组，按碰撞时间排序。
5. 在 mini_timestep 内发生的第一次碰撞会将您的 Universe（和 Universe 时钟）重置为 (last_time + time_to_collide)，其中 time_to_collide 没有更新。
6. 重复第 5 步，直到大 timestep 过去。按每个方格更新粒子的所有权。

该系统的优势在于，对于每个时间窗口，我们（假设粒子均匀分布）O(universe_particles * grid_size) 而不是 O(universe_particles * Universe_size) 检查碰撞。在良好的条件下（取决于宇宙大小、速度和粒子密度），您可以将计算效率提高几个数量级。

Answer 2

我不明白“优先队列”方法的工作原理，但我有一个替代方法可以帮助您。这就是我认为@Boyko Perfanov 的“利用本地化”的意思。

您可以将粒子分类到“桶”中，这样​​您就不必检查每个粒子之间的相互关系（ O(n²) ）。这利用了这样一个事实，即粒子只有在它们已经非常接近的情况下才能发生碰撞。创建代表小面积/体积的桶，并填充当前在桶的面积/体积中的所有粒子（ O(n) 最坏情况）。然后检查桶内的所有粒子与桶中的其他粒子（ O(m*(n/m)²) 平均情况，m = 桶数）。桶需要重叠才能工作，否则您还可以检查相邻桶中的粒子。

更新：如果粒子可以移动比桶大小更长的距离，一个明显的“解决方案”是减少时间步长。但是这会再次增加算法的运行时间，并且只有在有最大速度的情况下才有效。

即使没有最大速度也适用的另一种解决方案是创建一个额外的“高速”存储桶。由于速度分布通常是高斯曲线，因此不必将很多粒子放入该桶中，因此“桶方法”仍然比 O(n²) 更有效。