Hoomd-blue一款分子动力学模拟软件，由University of Michigan开发，开源于github，学习其详细内容可参考它的官方文档和github实例。

开端

题外话

Ethan在此之前计算机只有python基础课水平，没有独立参加过项目，所以学习了较长时间（断断续续有一个月），才将其基础搞懂，开始复现一篇论文中的结果，这篇论文（A geometry-originated universal relation for arbitrary convex hard particles）是导师要求学习的，但是没有想到看似简单的文章要学习和复现难度也如此大。frustrating

学习

论文主体上围绕一个公式展开

我也大概理解公式的意思，不过可能这也足够。左右侧都表征某种密度分布，左侧是用插入粒子来体现，右侧用多少粒子靠拢的程度来表征，公式的详细推导在论文里有。

复现准备

经过对官方文档孜孜不倦的学习，开始尝试写代码，（Hoomd文档没有类似的例子，这代码是真难写啊）。自认为代码非常粗糙，是Ethan在半知半解的调用函数和疯狂骚扰ChatGpt之后弄出来的结果：>

插入成功几率函数

为了得到插入成功的概率就需要一个能连续向体系中插入单个粒子的函数，不懈努力后，成功生造了一个在体系中插入粒子的函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def insert_triangle(simulation, mc,times=100):
    """
    在现有的 HOOMD 模拟中随机插入指定数量的三角形粒子，并检测与已有粒子是否重叠。

    参数：
    simulation: hoomd.Simulation 对象，表示现有的模拟。
    mc: hoomd的hpmc，蒙特卡洛积分器
    times：重复插入的次数，次数越大准确度越高
    """
    # 定义三角形的顶点，中心在 (0, 0)
    #a = 1.0  # 边长
    #h = np.sqrt(3) / 2 * a  # 等边三角形的高
    #vertices = [(-a/2, -h/3),
                #(a/2, -h/3),
                #(0, 2*h/3)]

    # 设置 HPMC 积分器
    #mc = hoomd.hpmc.integrate.ConvexPolygon()
    #mc.shape['A'] = {'vertices': vertices}

    new_mc = hoomd.hpmc.integrate.ConvexPolygon()
    new_mc.shape["A"] = dict(
        vertices = [
        (-0.5, -0.5),
        (0.5, -0.5),
        (0, 0.5),
    ]
    )

    # 将积分器添加到模拟中
    simulation.operations.integrator = mc
    

    # 获取模拟盒尺寸
    box = simulation.state.box
    Lx = box.Lx
    Ly = box.Ly

    # 获取粒子类型的索引
    type_id = simulation.state.particle_types.index('A')

    # 开始插入粒子
    inserted_recorder=0
    attempts = 0
    old_snap = simulation.state.get_snapshot()

    while attempts < times:
        attempts += 1
        # 随机生成位置
        x = np.random.uniform(-Lx/2, Lx/2)
        y = np.random.uniform(-Ly/2, Ly/2)
        #print(x,y)

        # 随机生成取向（四元数表示）
        #theta = np.random.uniform(0, 2*np.pi)
        #orientation = [np.cos(theta/2), 0, 0, np.sin(theta/2)]

        # 插入粒子
        if old_snap.communicator.rank == 0:
            # 记录旧的粒子数量
            N_old = old_snap.particles.N
            # 新的粒子数量
            N_new = N_old + 1
            theta = np.random.uniform(0, 2*np.pi)
            
            # 创建一个新的快照，具有更多的粒子
            new_snap = hoomd.Snapshot()
            new_snap.particles.N = N_new
            
            # 复制盒子尺寸
            new_snap.configuration.box = old_snap.configuration.box
            
            # 初始化粒子类型
            new_snap.particles.types = old_snap.particles.types
            
            # 初始化属性数组
            new_snap.particles.position[:] = np.zeros((N_new, 3), dtype=float)
            new_snap.particles.orientation[:] = np.zeros((N_new, 4), dtype=float)
            new_snap.particles.typeid[:N_old] = old_snap.particles.typeid[:]
            
            # 复制旧的粒子数据
            new_snap.particles.position[:N_old] = old_snap.particles.position[:]
            new_snap.particles.orientation[:N_old] = old_snap.particles.orientation[:]
            new_snap.particles.typeid[:N_old] = old_snap.particles.typeid[:]
            
            # 设置新粒子的属性
            new_snap.particles.position[N_old] = [x, y, 0]
            new_snap.particles.orientation[N_old] = [np.cos(theta/2), 0, 0, np.sin(theta/2)]
            new_snap.particles.typeid[N_old] = 0

            
            # 如果有其他属性（如 charge、diameter 等），也需要进行同样的处理
            
        else:
            new_snap = None  # 非主进程设置为 None

        cpu = hoomd.device.CPU()
        new_simulation = hoomd.Simulation(device=cpu,seed=1)
        new_simulation.create_state_from_snapshot(new_snap)
        new_simulation.operations.integrator = new_mc
        #check_snapshot=simulation.state.get_snapshot()
        #render(check_snapshot)
        #print(new_snap.particles.typeid)
        #print(new_snap.particles.position)

        # 检查重叠
        new_simulation.run(0)
        if new_mc.overlaps > 0:
            #print(f"检测到重叠，移除粒子 {inserted}")
            #remove_particle(simulation, simulation.state.N_particles - 1)
            #print(new_mc.overlaps)
            #print(new_mc.shape)
           
            #print(mc.overlaps)
            continue
        else:
            inserted_recorder += 1
        simulation.state.set_snapshot(old_snap)
        

    return inserted_recorder

需要了解的是，在hoomd中一般有两种储存粒子状态的方式，一种是GSD文件，方便更改，可以从外部轻松访问;另一种是snapshot，多用于在代码环境中储存系统信息，方便调用和删除释放内存。

SDF函数

SDF是scale distribution function的简称，中文尺度分布函数，毫无头绪。它是怎么一回事呢，假如体系粒子间的间隔同等乘以一个系数$1-x$（粒子之间的距离是中心与中心的距离）,那么有一定数量的粒子相碰，就可以计正好相碰（只有正好在相碰的粒子才算数）的粒子数为，好在这个函数是hoomd官方有的函数，不用现成写，只需要调用就行了

1
2
3
sdf_compute = hoomd.hpmc.compute.SDF(xmax=xmax, dx=dx)
sdf_xcompression = sdf_compute.x_compression
sdf_compression = sdf_compute.sdf_compression

此处SDF有两个重要的量可以通过类访问，x_compression和sdf_compression，它们分别是收缩时每一步的x和该x对应的sdf（x），两者都是numpy的array，所以我们后续很方便的就可以进行绘图

系统初始化

这一部分较简单且在hoomd的文档中花大篇幅和多个例子进行说明了，这里就不赘述原理，只针对本论文要求给出初始化函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def generate_particle(N_particles):
    #生成一个有N_particles的有序粒子群
    spacing = 2
    K = math.ceil(N_particles ** (1 / 2))
    L = K * spacing
    x = np.linspace(-L / 2, L / 2, K, endpoint=False)
    position_2d = list(itertools.product(x, repeat=2))  # 生成二维网格上的粒子位置

    # 取前 N_particles 个位置
    positions_2d = np.array(position_2d[0:N_particles])  # 形状为 (N_particles, 2)

    # 添加第三个坐标（z 坐标），设为零
    z_coordinates = np.zeros((N_particles, 1))  # 形状为 (N_particles, 1)
    positions_3d = np.hstack((positions_2d, z_coordinates))  # 合并为 (N_particles, 3)

    # 创建 GSD 帧并设置粒子属性
    frame = gsd.hoomd.Frame()
    frame.particles.N = N_particles
    frame.particles.position = positions_3d  # 现在是 (N_particles, 3) 的数组
    frame.particles.typeid = [0] * N_particles
    frame.configuration.box = [L, L, 0, 0, 0, 0]  # 注意这里的盒子尺寸

    #创建一个快照
    snapshot=hoomd.Snapshot()
    snapshot.particles.N = N_particles
    snapshot.particles.position[:] = positions_3d
    snapshot.particles.typeid[:] = [0] * N_particles
    snapshot.particles.types=['A']
    snapshot.configuration.box=[L, L, 0, 0, 0, 0]


    return snapshot

生成一个具有N_particles个粒子的体系，密度恒定

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def randomizing_particles(snapshot,times=10e3):
    cpu = hoomd.device.CPU()
    simulation = hoomd.Simulation(device=cpu, seed=29)

    mc = hoomd.hpmc.integrate.ConvexPolygon(default_d=0.3,default_a=0.4)
    mc.shape["A"] = dict(
        vertices = [
        (-0.5, -0.5),
        (0.5, -0.5),
        (0, 0.5),
    ]
    )

    simulation.operations.integrator = mc
    simulation.create_state_from_snapshot(snapshot)
    initial_snapshot = simulation.state.get_snapshot()

    simulation.run(times)
    final_snapshot = simulation.state.get_snapshot()

    return final_snapshot,mc,simulation

随机化这个粒子体系，times是MC随机化的步数

顺便一提，hoomd的官方文档的图像生成函数是基于 fresnel的三维图像生成函数。功能很强，但是我要做二维的体系，用不上：<

我就自己基于matplotlib写了一个二维图像生成函数（目前只支持展现粒子的位置，都是圆形的粒子）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def render(snapshot, dims=2):
    """
    从gsd文件生成一个二维的粒子群图像,目前仅支持圆粒子图.
    """
    # Extract particle positions and convert to a NumPy array
    positions = np.asarray(snapshot.particles.position)

    # Check dimensions (only 2D is supported in this function)
    if dims != 2:
        print("Only 2D rendering is supported in this function.")
        return

    # Plot particles
    plt.figure(figsize=(8, 8))
    plt.scatter(positions[:, 0], positions[:, 1], s=10, alpha=0.6)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Particle positions from GSD snapshot')
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

模拟与调教

进行中！！