HOOMD-blue学习笔记

分子模拟

2026-06-20

HOOMD-blue介绍

HOOMD-blue是一个Python软件包，可用于在CPU和GPU上对粒子系统进行模拟。它能够对各种形状的粒子进行蒙特卡罗模拟，同时也能对具有不同相互作用势的粒子进行分子动力学模拟。该软件的许多功能都是为软物质研究领域设计的，不过其代码本身具有通用性，可以用于多种类型的粒子模拟。

安装

mamba install conda-forge::hoomd

Tutorials 教程

1. 仿真对象

HOOMD-blue是一个面向对象的Python软件包。首先，需要导入该软件包：

import hoomd

模拟对象汇集了模拟所涉及的所有元素，并提供了用于运行模拟的接口。它由模拟状态以及作用于该状态的各类操作组成。模拟状态包括当前的区域、化学键、粒子的位置、速度、方向以及其他粒子属性。而各种操作则用于查看或修改这些状态。一个模拟只有一个状态，但可以包含任意数量的操作。

在构建模拟时，必须指定一个设备。该设备决定了模拟状态应存储在何处，以及执行操作时应使用哪款处理器。HOOMD-blue可以在CPU上运行：

cpu = hoomd.device.CPU()

或者GPU：

gpu = hoomd.device.GPU()

现在，你可以使用选定的设备来创建一个仿真实例。

simulation = hoomd.Simulation(device=cpu)

2. 硬粒子蒙特卡洛模拟

样板代码（Boilerplate Code）

import hoomd

粒子形状（Particle Shape）

硬粒子蒙特卡洛（HPMC）模拟将粒子表示为不可重叠的扩展物体。系统中不存在吸引或排斥力，粒子之间的相互作用完全由其几何形状决定。形式上讲，当粒子之间没有重叠时，系统的势能为零；一旦发生重叠，势能则为无穷大。这种纯粹的硬核相互作用会诱导粒子之间产生有效的“吸引力”，从而可能导致有序结构的形成。例如，规则的硬八面体会自组装成体心立方（bcc）结构。在本教程中，你将学习如何运行硬八面体的模拟，并观察这一行为。

八面体的自组装

积分器（The Integrator）

ConvexPolyhedron 积分器用于实现 HPMC 模拟。创建该积分器如下：

mc = hoomd.hpmc.integrate.ConvexPolyhedron()

通过设置 shape 属性来定义粒子的形状。一个凸多面体由一组顶点的凸包（convex hull）定义：

mc.shape["octahedron"] = dict(
    vertices=[
        (-0.5, 0, 0),
        (0.5, 0, 0),
        (0, -0.5, 0),
        (0, 0.5, 0),
        (0, 0, -0.5),
        (0, 0, 0.5),
    ]
)

试探移动（Trial Moves）

在每个时间步中，HPMC 会对系统中的每个粒子尝试 nselect 次试探移动。每次试探移动从伪随机数流中抽取，可能是平移移动（translation move）或旋转移动（rotation move）：

平移移动：将粒子沿随机方向移动一段随机距离（最大为 d）。
旋转移动：将粒子绕随机轴旋转一个随机角度。较大的 a 值允许更大的旋转幅度。

任何导致粒子形状与其他粒子重叠的试探移动都会被拒绝，粒子的位置和朝向保持不变；而不会引起重叠的试探移动则会被接受，粒子的位置或朝向将更新为新值。

nselect、d 和 a 都是积分器的属性，可按如下方式设置：

mc.nselect = 2
mc.d["octahedron"] = 0.15
mc.a["octahedron"] = 0.2

设置模拟（Setting the Integrator）

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)

积分器是一种操作（operation）。一个 Simulation 中只能有一个积分器，并且它会在每个时间步对系统状态进行操作。将 HPMC 积分器分配给模拟对象即可使用：

simulation.operations.integrator = mc

上面传入 Simulation 构造函数的 seed 参数用于选择伪随机数流中的数值序列。只要初始条件和种子相同，HPMC 模拟将产生完全相同的结果。

所有生成伪随机数的操作都使用模拟中设置的 seed。因此，每当添加使用随机数的操作时，都应将 seed 设置为非默认值。

现在，你已经创建了一个 Simulation 对象和一个 ConvexPolyhedron 积分器，但还不能运行模拟——你首先需要定义积分器所操作的系统状态。本教程的下一节将展示如何初始化系统状态。

3. 初始化系综状态

样板代码（Boilerplate Code）

import itertools
import math

import hoomd
import numpy

下一节（隐藏单元格）中的 render 函数使用 fresnel 渲染初始构型。源代码可在 hoomd-examples 仓库中找到。

系统状态的组成部分（Components of the System State）

在运行模拟之前，你需要先初始化系统状态。初始条件描述了模拟开始时系统中每个粒子的位置和朝向，以及周期性模拟盒子的尺寸。

本教程中的硬规则八面体系统会自组装成体心立方（bcc）结构。自组装是一种粒子在平衡状态下自发形成有序结构的过程。大多数自组装研究需要对数千个粒子进行数十小时的模拟。为了使本教程简短，我们仅模拟与 bcc 结构相容的少量粒子（数量为 $2 \times m^3$ ，其中 $m$ 为整数）。

m = 4
N_particles = 2 * m**3

放置粒子（Placing Particles）

在硬粒子蒙特卡洛（HPMC）中，有效的粒子构型不能存在重叠。本教程中的八面体粒子可内接于直径为 1 的球体内，因此我们将粒子放置在边长为 $L$ 的 $K \times K \times K$ 简单立方晶格上，间距略大于 1。

spacing = 1.2
K = math.ceil(N_particles ** (1 / 3))
L = K * spacing

在 HOOMD 中，粒子位置必须位于周期性盒子内部。立方盒子的范围是从 $-L/2$ 到 $L/2$ 。使用 itertools 和 numpy 将粒子位置放置在此范围内的晶格点上：

x = numpy.linspace(-L / 2, L / 2, K, endpoint=False)
position = list(itertools.product(x, repeat=3))
print(position[0:4])
# 输出示例：
# [(-3.6, -3.6, -3.6), (-3.6, -3.6, -2.4), (-3.6, -3.6, -1.2), (-3.6, -3.6, 0.0)]

由于 $K^3 \geq N_{\text{particles}}$ ，需将列表截断为恰好 $N_{\text{particles}}$ 个粒子：

position = position[0:N_particles]

你还需要为每个粒子设置朝向。HOOMD 使用四元数（quaternion）表示朝向。四元数 (1, 0, 0, 0) 表示无旋转。这里将所有粒子的朝向设为该值：

orientation = [(1, 0, 0, 0)] * N_particles

此时系统的构型如下图所示：

单位制（Units）

HOOMD-blue 并不采用任何特定的真实物理单位制，而是使用一套内部自洽的单位系统，并可兼容多种单位制。例如，如果你选择长度单位为米（meter）、能量单位为焦耳（Joule）、质量单位为千克（kilogram），那么力的单位自然就是牛顿（Newton），速度单位就是米/秒（m/s）。在纳米尺度系统中，常用单位包括纳米（nm）、千焦/摩尔（kJ/mol）和原子质量单位（amu）。

在 HPMC 中，主要的单位是长度。虽然质量存在，但它在配分函数中被约去，不会影响模拟过程。在无热（athermal）HPMC 系统中，能量尺度本身并不重要——因为重叠构型的能量为无穷大，而有效构型的势能为零。然而，在导出单位中仍会隐式出现能量，例如压强：

[\text{pressure}] = \frac{[\text{energy}]}{[\text{length}]^3}

在 HPMC 中，默认 $kT = 1$ （能量单位）。你可以通过以下公式将 HPMC 中的压强 $P$ 转换为无量纲形式：

P^* = \frac{P \sigma^3}{kT}

其中 $\sigma$ 是系统中的代表性长度，例如粒子的直径或边长。

将构型写入文件系统（Writing the Configuration to the File System）

GSD 文件用于存储周期性盒子、粒子位置、朝向以及其他状态属性。使用 GSD 的 Python 包来写入该文件：

import gsd.hoomd

Frame 对象用于存储系统状态：

frame = gsd.hoomd.Frame()
frame.particles.N = N_particles
frame.particles.position = position
frame.particles.orientation = orientation

每个粒子还有一个类型（type）。在 HPMC 中，每种类型都有其独立的形状参数。本教程只有一种类型，因此将所有粒子的类型 ID 设为 0（类型 ID 从 0 开始编号）：

frame.particles.typeid = [0] * N_particles

每种粒子类型都需要一个名称。名称可以是任意字符串。HOOMD-blue 使用类型名称将你在操作中指定的参数与状态中的类型名称进行匹配。这里将类型 ID 0 命名为 'octahedron'：

frame.particles.types = ["octahedron"]

GSD 使用一个包含 6 个元素的数组表示盒子：三个边长 $L_x, L_y, L_z$ 和三个倾斜因子（tilt factors）。设置相等的边长和零倾斜因子以定义一个立方盒子：

frame.configuration.box = [L, L, L, 0, 0, 0]

将此快照写入 lattice.gsd 文件：

with gsd.hoomd.open(name="lattice.gsd", mode="x") as f:
    f.append(frame)

初始化 Simulation（Initializing a Simulation）

你可以使用该文件来初始化 Simulation 的状态。首先创建一个 Simulation 对象：

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu)

然后使用 create_state_from_gsd 工厂函数读取 GSD 文件，并用其中的粒子和盒子信息填充模拟状态：

simulation.create_state_from_gsd(filename="lattice.gsd")

本教程的下一节将展示如何使用 HPMC 对粒子的位置和朝向进行随机化处理。

4. 初始化系综状态

import math

import hoomd

方法（Method）

本教程的上一节将所有粒子放置在一个简单立方晶格上。这是一种方便的、能确保粒子不重叠的初始排布方式，但它使模拟从一个高度有序的状态开始。你应该对系统进行充分的随机化，使其“忘记”这一初始状态，从而让自组装过程不受初始构型的影响。

你不能简单地为粒子位置生成完全随机的数，因为那样会导致粒子之间发生重叠。取而代之的方法是：从晶格构型出发，利用 HPMC 在保证无重叠的前提下对粒子进行随机移动。在低密度构型下（如上一节生成的晶格），只需运行较短时间的模拟即可快速实现系统的随机化。

设置 Simulation（Set up the Simulation）

以下代码块创建了 Simulation 对象，配置了 HPMC 积分器，并从 lattice.gsd 文件中初始化了系统状态——这些内容已在本教程前面的章节中讨论过：

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=12)

mc = hoomd.hpmc.integrate.ConvexPolyhedron()
mc.shape["octahedron"] = dict(
    vertices=[
        (-0.5, 0, 0),
        (0.5, 0, 0),
        (0, -0.5, 0),
        (0, 0.5, 0),
        (0, 0, -0.5),
        (0, 0, 0.5),
    ]
)

simulation.operations.integrator = mc
simulation.create_state_from_gsd(filename="lattice.gsd")

运行模拟（Run the Simulation）

保存当前系统状态的一个快照。本教程稍后会用它来观察粒子移动了多远：

initial_snapshot = simulation.state.get_snapshot()

运行模拟以随机化粒子的位置和朝向。run 方法接受要运行的时间步数作为参数。对于低密度系统，10,000 步已足以完成充分的随机化：

simulation.run(10e3)

你可以查询 HPMC 积分器的属性，了解它执行了哪些操作：

translate_moves 是一个元组，包含被接受和平移被拒绝的移动次数。在此低密度下，平移移动的接受率（即被接受的试探移动占总试探移动的比例）非常高：

mc.translate_moves
# (2412396, 148877)

mc.translate_moves[0] / sum(mc.translate_moves)
# 0.9418738260232314

rotate_moves 同样提供了被接受和被拒绝的旋转移动次数：

mc.rotate_moves
# (2486740, 71987)

mc.rotate_moves[0] / sum(mc.rotate_moves)
# 0.9718660880977142

overlaps 报告当前状态中重叠粒子对的数量。最终构型中没有重叠：

mc.overlaps
# 0

最终构型（The Final Configuration）

查看最终的粒子位置和朝向，看看它们发生了怎样的变化：

final_snapshot = simulation.state.get_snapshot()
render(final_snapshot)

tutorial_00-Introducing-HOOMD-blue_04-Randomizing-the-System_21_0.png

初始位置（前4个粒子）：

initial_snapshot.particles.position[0:4]
# array([[-3.5999999 , -3.5999999 , -3.5999999 ],
#        [-3.5999999 , -3.5999999 , -2.4000001 ],
#        [-3.5999999 , -3.5999999 , -1.20000005],
#        [-3.5999999 , -3.5999999 ,  0.        ]])

最终位置（前4个粒子）：

final_snapshot.particles.position[0:4]
# array([[-0.5163401 , -1.72834745,  1.70890652],
#        [-3.49591272, -0.29685889,  2.04439483],
#        [-2.24041986, -2.3811753 , -2.62242222],
#        [-2.52071304,  2.813749  , -1.10457279]])

初始朝向（前4个粒子）：

initial_snapshot.particles.orientation[0:4]
# array([[1., 0., 0., 0.],
#        [1., 0., 0., 0.],
#        [1., 0., 0., 0.],
#        [1., 0., 0., 0.]])

最终朝向（前4个粒子）：

final_snapshot.particles.orientation[0:4]
# array([[ 0.27468776, -0.67511118, -0.66949342, -0.14335298],
#        [-0.80225824, -0.12547562,  0.25790363,  0.52356786],
#        [ 0.43336577, -0.28480202, -0.02206664, -0.85474849],
#        [ 0.44913914, -0.30817778, -0.75143279,  0.37235636]])

可以看到，粒子的位置和朝向确实发生了显著变化，说明系统已被充分随机化。

将最终构型保存为 GSD 文件，供模拟的下一阶段使用：

hoomd.write.GSD.write(state=simulation.state, mode="xb", filename="random.gsd")

本教程的下一节将使用 random.gsd 文件，并将其压缩至更高密度。

5. 压缩体系

体积分数（Volume Fraction）

硬粒子系统中的自组装通常发生在体积分数高于 0.5 的条件下。体积分数定义为粒子所占体积与周期性模拟盒子体积的比值。

到目前为止，本教程已在一个体积分数非常低的盒子中对 $N$ 个八面体进行了随机化，并将结果保存在 random.gsd 文件中。现在，用该构型初始化一个 Simulation，并查看当前的体积分数：

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=12)
simulation.create_state_from_gsd(filename="random.gsd")

计算单个八面体的体积：

a = math.sqrt(2) / 2
V_particle = (1 / 3) * math.sqrt(2) * a**3

计算当前体积分数：

initial_volume_fraction = (
    simulation.state.N_particles * V_particle / simulation.state.box.volume
)
print(initial_volume_fraction)
# 0.05715592589579699

如你所见，当前体积分数非常低，必须大幅减小盒子体积，才能达到高于 0.5 的目标体积分数。

我们将使用 HPMC 在压缩系统的同时，使粒子移动到无重叠的构型中。首先为八面体模拟设置 HPMC 积分器：

mc = hoomd.hpmc.integrate.ConvexPolyhedron()
mc.shape["octahedron"] = dict(
    vertices=[
        (-0.5, 0, 0),
        (0.5, 0, 0),
        (0, -0.5, 0),
        (0, 0.5, 0),
        (0, 0, -0.5),
        (0, 0, 0.5),
    ]
)
simulation.operations.integrator = mc

QuickCompress 更新器（The QuickCompress Updater）

在 HOOMD-blue 中，Updater（更新器） 是一种会修改系统状态的操作。要使用 Updater，需先实例化对象、分配一个 Trigger（触发器），然后将其添加到 Simulation 中。当 Trigger 返回 True 时，Simulation 就会在对应的时间步应用该 Updater。Periodic 触发器每隔固定步数返回一次 True。

QuickCompress 是一个专为 HPMC 设计的 Updater，可快速将盒子压缩至目标体积。每次被触发时，它会按比例缩小盒子体积，同时允许粒子之间出现轻微重叠；随后等待平移和旋转试探移动消除这些重叠，再进行下一次压缩。此过程会暂时产生无效构型（即存在重叠），但比完全禁止重叠的方法快得多。

现在，我们计算目标体积分数为 0.57 时所需的最终盒子体积，并配置一个每 10 步触发一次的 QuickCompress：

initial_box = simulation.state.box
final_box = hoomd.Box.from_box(initial_box)
final_volume_fraction = 0.57
final_box.volume = simulation.state.N_particles * V_particle / final_volume_fraction

compress = hoomd.hpmc.update.QuickCompress(
    trigger=hoomd.trigger.Periodic(10), target_box=final_box
)

将该 Updater 添加到 Simulation 中：

simulation.operations.updaters.append(compress)

MoveSize 调优器（The MoveSize Tuner）

Tuner（调优器） 是另一种操作类型，用于动态调整其他操作的参数以提升性能。在 HPMC 中，平移和旋转试探移动的步长（d 和 a）对性能影响极大：

若步长太小，系统需要很多时间步才能发生显著变化；
若步长太大，则绝大多数移动会被拒绝，同样导致效率低下。

通常，当接受率（acceptance ratio）约为 20% 时，系统演化效率最高。MoveSize 调优器会监控接受率，并自动调整 d（平移步长）和 a（旋转步长）以达到目标接受率。

由于系统密度在 QuickCompress 压缩过程中快速变化，最优步长也随之改变。因此，我们定期使用 MoveSize 调优器进行调整：

periodic = hoomd.trigger.Periodic(10)
tune = hoomd.hpmc.tune.MoveSize.scale_solver(
    moves=["a", "d"],
    target=0.2,
    trigger=periodic,
    max_translation_move=0.2,
    max_rotation_move=0.2,
)
simulation.operations.tuners.append(tune)

运行直至完成（Run Until Complete）

当 QuickCompress 达到目标盒子尺寸且系统中无任何粒子重叠时，压缩过程即告完成。所需时间步数取决于具体参数。我们通过循环定期检查 compress.complete 属性，并在压缩完成后停止模拟：

while not compress.complete and simulation.timestep < 1e6:
    simulation.run(1000)

其中 simulation.timestep < 1e6 的判断是为了防止在压缩无法完成的情况下无限循环浪费计算资源。正常情况下，循环会在远小于 100 万步时结束：

simulation.timestep
# 13000

检查压缩是否成功完成：

if not compress.complete:
    message = "Compression failed to complete"
    raise RuntimeError(message)

此时，MoveSize 调优器应已根据高密度环境将移动步长调整为较小的值：

mc.a["octahedron"]
# 0.0579430568204312

mc.d["octahedron"]
# 0.030998272372613312

现在压缩已完成，粒子彼此非常接近但无重叠，且仍保持随机排布。

render(simulation.state.get_snapshot())

tutorial_00-Introducing-HOOMD-blue_05-Compressing-the-System_30_0.png

最后，将最终构型保存为 GSD 文件，供模拟的下一阶段使用：

hoomd.write.GSD.write(state=simulation.state, mode="xb", filename="compressed.gsd")

6. 平衡系统

以下是您提供英文内容的中文翻译：

平衡化（Equilibration）

到目前为止，本教程已经将 $N$ 个不重叠的八面体随机放置在一个盒子中，并将其压缩至中等体积分数。所得的粒子构型虽然是有效的，但强烈依赖于生成它的路径。而在所有可能的构型集合中，存在大量与路径无关的平衡态构型。平衡化（Equilibration） 的过程就是从一个人为制备的初始状态出发，运行模拟，使系统在演化过程中逐渐弛豫至热力学平衡态。

首先初始化 Simulation：

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=12)

mc = hoomd.hpmc.integrate.ConvexPolyhedron()
mc.shape["octahedron"] = dict(
    vertices=[
        (-0.5, 0, 0),
        (0.5, 0, 0),
        (0, -0.5, 0),
        (0, 0.5, 0),
        (0, 0, -0.5),
        (0, 0, 0.5),
    ]
)
simulation.operations.integrator = mc

本教程上一节已将压缩后的系统保存为 compressed.gsd。现在从该文件初始化系统状态：

simulation.create_state_from_gsd(filename="compressed.gsd")

写入模拟轨迹（Writing Simulation Trajectories）

为了观察平衡化过程，需定期将系统状态保存到文件中。此前本教程使用 GSD 文件通过 GSD Python 包或 GSD.write 保存单帧状态。现在我们将使用 GSD Writer（另一种操作） 来创建一个包含多帧的轨迹文件（trajectory）。

我们使用 'xb' 模式打开文件，以确保在写入前该文件不存在。如果文件已存在，则会抛出错误，而不是覆盖或追加内容：

gsd_writer = hoomd.write.GSD(
    filename="trajectory.gsd", trigger=hoomd.trigger.Periodic(1000), mode="xb"
)
simulation.operations.writers.append(gsd_writer)

调优试探移动步长（Tuning the Trial Move Size）

上一节在压缩过程中频繁使用了 MoveSize 调优器，在系统密度快速变化时动态调整 d（平移步长）和 a（旋转步长），以维持目标接受率。在平衡化阶段，我们也应再次使用它，以确保 HPMC 高效运行。

不过需要注意：移动步长应在模拟初期短暂调优，之后在整个运行过程中保持恒定。若在整个模拟过程中持续改变步长，会破坏细致平衡（detailed balance），导致结果不正确。

为此，我们将调优器设置为：仅在前 5000 步内，每隔 100 步触发一次。这通过组合 Periodic 和 Before 触发器，并用 And 逻辑实现：

Before(t) 在时间步 $t < \text{value}$ 时返回 True；
And 触发器仅在其所有子触发器都返回 True 时才返回 True。

tune = hoomd.hpmc.tune.MoveSize.scale_solver(
    moves=["a", "d"],
    target=0.2,
    trigger=hoomd.trigger.And(
        [hoomd.trigger.Periodic(100), hoomd.trigger.Before(simulation.timestep + 5000)]
    ),
)
simulation.operations.tuners.append(tune)

运行前 5000 步进行调优：

simulation.run(5000)

接着再运行 100 步，检查接受率是否接近目标值（约 20%）：

simulation.run(100)

rotate_moves = mc.rotate_moves
mc.rotate_moves[0] / sum(mc.rotate_moves)
# 0.20826168516492105

translate_moves = mc.translate_moves
mc.translate_moves[0] / sum(mc.translate_moves)
# 0.19951585194440105

可以看到，平移和旋转的接受率均已非常接近 20%，说明调优成功。

系统平衡化（Equilibrating the System）

要使系统达到平衡，只需继续运行模拟。所需的时间步数高度依赖于具体模型、系统规模、密度等多种因素。对于约 10,000 个粒子的硬粒子蒙特卡洛自组装系统，通常需要数千万步才能充分平衡。而本教程中的系统规模小得多，因此所需步数也更少。

⏱️ 注意：以下单元格可能需要几分钟才能完成。

simulation.run(2e5)  # 运行 200,000 步

这是运行结束后系统的最终状态：

render(simulation.state.get_snapshot())

tutorial_00-Introducing-HOOMD-blue_06-Equilibrating-the-System_21_0.png

这个最终状态是否为平衡态？ 本教程的下一节将展示如何分析轨迹文件（trajectory），并回答这一问题。

关闭 GSD 文件

如果你在不关闭本 notebook 的情况下直接运行下一节的内容，必须手动关闭 GSD 文件。通常在 Python 脚本结束时 GSD 文件会自动关闭，但在 Jupyter 环境中需要显式处理：

simulation.operations.writers.remove(gsd_writer)
del gsd_writer

#分子模拟 ⁵

Haddock3学习笔记

小康文章阅读笔记