HOOMD-blue——Modelling Rigid Bodies

风往北吹

教程介绍了如何运行刚体的模拟。
刚体简介

1.导入包

import itertools
import math

import gsd.hoomd
import hoomd
import matplotlib
import numpy

%matplotlib inline
matplotlib.style.use('ggplot')
2.简介

刚体是一种不可压缩体，由一个中心粒子和一个或多个组分粒子组成。在MD模拟中，中心粒子根据作用于自身的力和扭矩进行平移和旋转，组分粒子跟随中心粒子运动。
3.坐标

我们在体坐标（body coordinates）中定义组分粒子的位置和取向，其中(0,0,0)是中心粒子的位置。例如，我们定义一个刚性二聚体中两个点的位置：
dimer_positions = [[-1.2, 0, 0], [1.2, 0, 0]]

每个刚体会根据中心粒子的位置和取向，被放置在全局坐标（global coordinates）。
central_position = [10, 5, 0]
central_rotation = 0.9

下面的代码计算了在给定中心位置和取向的二聚体在全局坐标中的位置（2D）：
cos_theta = math.cos(central_rotation)
sin_theta = math.sin(central_rotation)

global_positions = []
for i in range(len(dimer_positions)):
    x, y = dimer_positions[:2]

    global_positions.append(
        [[central_position[0] + (x * cos_theta - y * sin_theta)],
         [central_position[1] + (y * cos_theta + x * sin_theta)]])

我们可以进行可视化查看：
fig = matplotlib.figure.Figure(figsize=(10, 6.18), dpi=100)
ax = fig.add_subplot(aspect='equal')

ax.plot([central_position[0], central_position[0] - 3 * sin_theta],
        [central_position[1], central_position[1] + 3 * cos_theta],
        color='k')

ax.plot([central_position[0], central_position[0] + 3 * cos_theta],
        [central_position[1], central_position[1] + 3 * sin_theta],
        color='k')

ax.text(central_position[0] + 1.5 * cos_theta,
        central_position[1] + 1.5 * sin_theta,
        'Body X',
        rotation=central_rotation * 180 / math.pi,
        verticalalignment='center')

ax.text(central_position[0] - 2 * sin_theta,
        central_position[1] + 2 * cos_theta,
        'Body Y',
        rotation=central_rotation * 180 / math.pi + 90,
        verticalalignment='center')

ax.add_patch(
    matplotlib.patches.Circle((central_position[0], central_position[1]),
                              0.1,
                              color='C0'))

for position in global_positions:
    ax.add_patch(
        matplotlib.patches.Circle((position[0], position[1]), 1.0, color='C1'))

ax.set_xlim(0, 20)
ax.set_ylim(0, 12)

ax.set_xlabel('Global X')
ax.set_ylabel('Global Y')

fig


4.刚体的属性

刚体的属性主要表现在中心粒子的属性，主要包括：质量-刚体的总质量；惯性矩张量-质量如何分布在整个刚体上；速度-刚体质心的速度；角动量-刚体的角动量；位置-全局坐标中的位置；取向-四元数，使刚体围绕中心粒子旋转。

还包括组分粒子的属性：组分粒子的位置-体坐标中刚体中心到组分粒子的向量；组分粒子取向-组分粒子绕自身中心旋转的取向。
5.定义刚体二聚体的属性

我们将二聚体中的组分粒子设置为质量为1，类型为A的点粒子，位置为前述的dimer_positions。我们将每个二聚体的类型定义为dimer。

现在开始创建二聚体的初始构型。首先创建包含两种粒子类型的snapshot：
snapshot = gsd.hoomd.Snapshot()
snapshot.particles.types = ['dimer', 'A']

将二聚体（typeid=0）放置在模拟盒子中，将二聚体粒子之间的距离放置的足够远以便A粒子之间不会接触：
m = 4
N_particles = m**3
spacing = 5
K = math.ceil(N_particles**(1 / 3))
L = K * spacing
x = numpy.linspace(-L / 2, L / 2, K, endpoint=False)
position = list(itertools.product(x, repeat=3))
position = numpy.array(position) + [spacing / 2, spacing / 2, spacing / 2]
snapshot.particles.N = N_particles
snapshot.particles.position = position[0:N_particles, :]
snapshot.particles.typeid = [0] * N_particles
snapshot.configuration.box = [L, L, L, 0, 0, 0]

由于两个A粒子的质量均为1，所以一个二聚体的质量为2。设置二聚体的质量：
snapshot.particles.mass = [2] * N_particles

特定轴的点粒子的惯性矩为,为点离轴的距离。通常情况下，惯性矩为一个张量并且包含非对角线的值。

我们可以计算得到二聚体的惯性张量：
mass = 1
I = numpy.zeros(shape=(3, 3))
for r in dimer_positions:
    I += mass * (numpy.dot(r, r) * numpy.identity(3) - numpy.outer(r, r))
I
array([[0.  , 0.  , 0.  ],
       [0.  , 2.88, 0.  ],
       [0.  , 0.  , 2.88]])

在本例中，该张量是对角化的。在HOOMD-blue中假设刚体在体坐标中具有对角化的惯性矩，其形式为。我们可以在snapshots中设置每个刚体的惯性矩：
snapshot.particles.moment_inertia = [I[0, 0], I[1, 1], I[2, 2]] * N_particles

本例中是0，而和都是非零值。HOOMD检查不为零的惯性矩并对这些惯性矩的轴进行自由度的积分。在本例子中，二聚体将绕着刚体的y和z轴旋转，而不绕x轴进行旋转。

HOOMD用四元数表示取向。我们将每个刚体的四元数设置为1：
snapshot.particles.orientation = [(1, 0, 0, 0)] * N_particles

然后将刚体的中心写入到GSD文件中：
with gsd.hoomd.open(name='dimer_centers.gsd', mode='xb') as f:
    f.append(snapshot)

使用hoomd.md.contrain.Rigid对刚体进行约束：
rigid = hoomd.md.constrain.Rigid()

然后将组分粒子的属性设置给刚体：
rigid.body['dimer'] = {
    "constituent_types": ['A', 'A'],
    "positions": dimer_positions,
    "orientations": [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)],
    "charges": [0.0, 0.0],
    "diameters": [1.0, 1.0]
}
6.在初始状态中放置组分粒子

到目前位置，我们的snapshot中只有二聚体的中心粒子位置：
sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=4)
sim.create_state_from_gsd(filename='dimer_centers.gsd')
render(sim.state.get_snapshot())


Rigid.create_bodies能够将组分粒子放置在盒子中：
rigid.create_bodies(sim.state)

现在二聚体的每个中心粒子都有两个组分粒子：

要在MD积分器中设置rigid约束才能使其在模拟中起作用：
integrator = hoomd.md.Integrator(dt=0.005, integrate_rotational_dof=True)
integrator.rigid = rigid
sim.operations.integrator = integrator

Rigid将在模拟开始后的每一步对刚体进行约束。调用sim.run(0)将储存组分粒子的位置：
sim.run(0)
render(sim.state.get_snapshot())


要想改变刚体的位置和取向，需要对中心粒子的属性进行改变：
with sim.state.cpu_local_snapshot as snapshot:
    typeid = snapshot.particles.typeid
    snapshot.particles.orientation[typeid == 0] = [
        0.70710678, 0., 0.70710678, 0.
    ]

然后再执行sim.run(0)使改变生效：
sim.run(0)
render(sim.state.get_snapshot())


然后将初始构型写入到GSD文件中以便后续使用：
hoomd.write.GSD.write(state=sim.state, mode='xb', filename='lattice.gsd')
运行刚体模拟

1.导入包

import math

import hoomd
2.从GSD文件继续进行模拟

上述部分我们创建了包含刚体二聚体初始构型的lattice.gsd文件。现在我们从这个文件中初始化一个新的模拟：
cpu = hoomd.device.CPU()
sim = hoomd.Simulation(device=cpu, seed=1)
sim.create_state_from_gsd(filename='lattice.gsd')

然后我们还需要构建刚体约束：
rigid = hoomd.md.constrain.Rigid()

刚体的体属性（组分粒子的属性）并没有保存在GSD文件中，因此我们需要重新定义一遍：
rigid.body['dimer'] = {
    "constituent_types": ['A', 'A'],
    "positions": [[-1.2, 0, 0], [1.2, 0, 0]],
    "orientations": [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)],
    "charges": [0.0, 0.0],
    "diameters": [1.0, 1.0]
}
3.积分刚体的自由度

使用integrate_rotational_dof=True选项构建积分器：
integrator = hoomd.md.Integrator(dt=0.005, integrate_rotational_dof=True)
sim.operations.integrator = integrator

默认情况下integrate_rotational_dof为False。当值为False时，积分器将会忽略所有粒子的扭力，只对平移自由度进行积分。当为True时，积分器对平移和旋转自由度均进行积分。

将刚体约束分配给积分器属性：
integrator.rigid = rigid

平移和旋转自由度是赋给刚体的中心粒子。在积分过程中，只选择刚体中心和自由粒子进行积分：
kT = 1.5
rigid_centers_and_free = hoomd.filter.Rigid(("center", "free"))
nvt = hoomd.md.methods.NVT(kT=kT, tau=1., filter=rigid_centers_and_free)
integrator.methods.append(nvt)
4.定义刚体之间的对相互作用

我们需要创建一个带有body exclusion的邻居列表，这样将不会计算刚体内部的相互作用：
cell = hoomd.md.nlist.Cell(buffer=0, exclusions=['body'])

如果没有body exclusion，模拟可能会因为大数相减时的四舍五入误差而在数值上不稳定。不管在哪种情况，刚体内部的相互作用均对力和扭力没有贡献，因此不需要计算他们。

在A-A粒子之间使用LJ势：
lj = hoomd.md.pair.LJ(nlist=cell)
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)
lj.r_cut[('A', 'A')] = 2.5

该体系中还存在中心粒子（type dimer），只要将dimer粒子设置为无相互作用即可：
lj.params[('dimer', ['dimer', 'A'])] = dict(epsilon=0, sigma=0)
lj.r_cut[('dimer', ['dimer', 'A'])] = 0

然后将该势能加入到积分器中：
integrator.forces.append(lj)
5.压缩体系

我们已经对体系状态进行了初始化，添加了积分器，设置了刚体约束并定义了粒子对相互作用。在运行模拟之前，还需要设置初始速度和刚体中心粒子和自由粒子的角动量：
sim.state.thermalize_particle_momenta(filter=rigid_centers_and_free, kT=kT)
sim.run(0)
nvt.thermalize_thermostat_dof()

thermalize_particle_momenta给filter中的每个粒子的每个平移和旋转自由度分配高斯分布的动量。在本例中，将会给y和z轴分配动量，而不给x轴分配因为。

现在我们可以运行随机化的模拟：
sim.run(1000)
render(sim.state.get_snapshot())


然后对体系进行压缩：
initial_box = sim.state.box
final_box = hoomd.Box.from_box(initial_box)  # make a copy of initial_box
final_box.volume = initial_box.volume / 16
box_resize_trigger = hoomd.trigger.Periodic(10)
ramp = hoomd.variant.Ramp(A=0, B=1, t_start=sim.timestep, t_ramp=20000)
box_resize = hoomd.update.BoxResize(box1=initial_box,
                                    box2=final_box,
                                    variant=ramp,
                                    trigger=box_resize_trigger)
sim.operations.updaters.append(box_resize)

sim.run(30000)
render(sim.state.get_snapshot())


6.刚体的热力学性质

我们可以添加ThermodynamicQuantities，在模拟中计算热力学性质：
thermodynamic_quantities = hoomd.md.compute.ThermodynamicQuantities(
    filter=hoomd.filter.All())

sim.operations.computes.append(thermodynamic_quantities)

ThermodynamicQuantities分别计算了平移和旋转动能：
thermodynamic_quantities.translational_kinetic_energy
136.95366173865793
thermodynamic_quantities.rotational_kinetic_energy
105.62170571721094
准备一个任意的刚体

1.导入包

import math

import numpy
2.刚体惯性张量对角化

HOOMD-blue假设体系中提供的刚体的惯性矩张量是对角化的。

给定一个任意放置组分粒子的刚体，可以使用下面的方法来准备可在HOOMD中使用的刚体形式。

第一步，得到刚体的惯性张量（非对角化的）；

第二步，求解使惯性张量对角化的刚体的取向；

第三步，将刚体重新取向使得其惯性张量在全局坐标中使对角化的；

第四步，使用新的对角化的惯性矩。如下面的例子所示，首先我们定义任意的刚体：
general_positions = numpy.array([[.5, .5, 0], [-.5, -.5, 0], [-1, 1, 0],
                                 [1, -1, 0]])
particle_mass = 1
particle_radius = 1


1.计算惯性张量

HOOMD-blue对于刚体的质量分布不做任何假设。例如，假设每个组分粒子是一个密度均匀的球，我们可以计算惯性矩。单个球围绕其中心的惯性矩为：
I_ref = numpy.array([[2 / 5 * particle_mass * particle_radius**2, 0, 0],
                     [0, 2 / 5 * particle_mass * particle_radius**2, 0],
                     [0, 0, 2 / 5 * particle_mass * particle_radius**2]])
I_ref
array([[0.4, 0. , 0. ],
       [0. , 0.4, 0. ],
       [0. , 0. , 0.4]])

利用平行轴定理，计算一般物体的惯性张量为各组分粒子的惯性张量之和：。

是参考的惯性矩，m为质量，D为位移（即组分粒子在体坐标中的位置），为单位矩阵。
I_general = numpy.zeros(shape=(3, 3))
for r in general_positions:
    I_general += I_ref + particle_mass * (numpy.dot(r, r) * numpy.identity(3)
                                          - numpy.outer(r, r))

I_general

array([[4.1, 1.5, 0. ],
       [1.5, 4.1, 0. ],
       [0. , 0. , 6.6]])

我们可以发现一些非对角的元素是非零的。
2.惯性矩对角化

首先找到该惯性矩的特征值：
I_diagonal, E_vec = numpy.linalg.eig(I_general)

在HOOMD-blue中所需要的对角化的惯性矩为。
I_diagonal
array([5.6, 2.6, 6.6])

特征向量矩阵的转置能够将原始粒子旋转到惯性张量被对角化的框架中：
R = E_vec.T

使用该矩阵对粒子进行旋转：
diagonal_positions = numpy.dot(R, general_positions.T).T
diagonal_positions
array([[ 7.07106781e-01, -5.55111512e-17,  0.00000000e+00],
       [-7.07106781e-01,  5.55111512e-17,  0.00000000e+00],
       [ 1.11022302e-16,  1.41421356e+00,  0.00000000e+00],
       [-1.11022302e-16, -1.41421356e+00,  0.00000000e+00]])
render_rigid_body(diagonal_positions)



现在完成了对角化后，我们可以使用diagonal_positions作为体坐标，同时用I_diagonal作为中心粒子的惯性矩。

现在可以对惯性矩进行检查是否对角化成功：
I_check = numpy.zeros(shape=(3, 3))
for r in diagonal_positions:
    I_check += I_ref + particle_mass * (numpy.dot(r, r) * numpy.identity(3)
                                        - numpy.outer(r, r))

I_check
array([[ 5.60000000e+00, -2.35513869e-16,  0.00000000e+00],
       [-2.35513869e-16,  2.60000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  6.60000000e+00]])

详细教程可参考HOOMD-blue官方网站：

https://hoomd-blue.readthedocs.io/en/v3.6.0/tutorial/06-Modelling-Rigid-Bodies/00-index.html