周期边界盒中的粒子-粒子碰撞

Question

对于分子动力学的一门课，我要模拟一个有周期边界的盒子里的100个粒子。我必须考虑粒子与粒子的碰撞，因为墙壁是“透明的”，这些相互作用可以跨越边界发生。由于模拟应该包含50000个步骤，而且我希望将来会有越来越多的额外任务，所以我希望我的代码尽可能高效（尽管运行时间很长，但我必须使用python）。

该系统包括

null

class particlesInPeriodicBox:
    # define the particle properties
    def __init__(self,
                 initialState = [[0,0,0,0]], # state with form [x, y, vx, vy] for each particle
                 boundaries = [0, 10, 0, 5], # box boundaries with form [xmin, xmax, ymin, ymax] in nm
                 radius = 0.2, # particle radius in nm
                 mass = 2): # mass in g/mol. Here a parameter instead of a global variable
        self.initialState = np.asarray(initialState, dtype=float)
        self.state = self.initialState.copy()
        self.radius = radius
        self.time = 0 # keep count of time, if time, i want to implement a 'clock'
        self.mass = mass # mass
        self.boundaries = boundaries
    def collision(self):
        """
        now, one has to check for collisions. I do this by distance check and will solve this problems below.
        To minimize the simulation time I then will only consider the particles that colided.
        """
        dist = squareform(pdist(self.state[:, :2])) # direct distance

        colPart1, colPart2 = np.where(dist < 2 * self.radius) # define collision partners 1 and 2 as those where the distance between the centeres are smaller than two times the radius

        # resolve self-self collissions
        unique = (colPart1 < colPart2)
        colPart1 = colPart1[unique]
        colPart2 = colPart2[unique]

        """
        The following loop resolves the collisions. I zip the lists of collisionpartners to one aray,
        where one entry contains both colisionpartners.
        """
        for cp1, cp2 in zip(colPart1, colPart2): # cp1/cp2 are the two particles colliding in one collision.
            # masses could be different in future...
            m1 = self.mass[cp1]
            m2 = self.mass[cp2]

            # take the position (x,y) tuples for the two particles
            r1 = self.state[cp1, :2]
            r2 = self.state[cp2, :2]

            # same with velocities
            v1 = self.state[cp1, 2:]
            v2 = self.state[cp2, 2:]

            # get relative parameters
            r = r1-r2
            v = r2-r1

            # center of mass velocity:
            vcm = (m1 * v1 + m2 * v2) / (m1 + m2)

            """
            This is the part with the elastic collision
            """
            dotrr = np.dot(r, r) # the dot product of the relative position with itself
            dotvr = np.dot(v, r) # the dot product of the relative velocity with the relative position 
            v = 2 * r * dotvr / dotrr - v # new relative velocity

            """
            In this center of mass frame, the velocities 'reflect' on the center of mass in oposite directions
            """
            self.state[cp1, 2:] = vcm + v * m2/(m1 + m2) # new velocity of particle 1 still considering possible different masses
            self.state[cp2, 2:] = vcm - v * m1/(m1 + m2) # new velocity of particle 2

按照我的理解，这种处理对整个数组的操作的技术比每次手动遍历它更有效。移动的粒子‘槽’的墙壁是容易的，我只是减去或增加的盒子的尺寸，分别。但是：

目前，该算法只看到方框内部的冲突，而不是跨越边界。这个问题我想了一会，想出了以下几个想法：

• 我可以在一个3x3网格中制作该系统的总共9个副本，并且只考虑中间的一个，因此可以查看相邻单元格以进行最近邻搜索。但是我想不出一个有效的方法来实现它，尽管这个方法是标准方法
• 其他的想法都有模数的使用，我几乎可以肯定，这不是正确的方法。

null

对于第一个问题，也许可以使用一些技术，如在周期性边界条件下计算接触数/配位数，但我不确定这是否是最有效的方法。

谢谢你！

Answer 1

模数很可能是一个快速的运算，就像你要得到的一样。在任何自重的运行时系统中，这将直接连接到片内浮除操作，这可能比一组繁琐的“if-subtract”对更快。

我觉得你的9细胞解决方案有点过火了。在2x2矩阵中使用4个单元格，检查两个区域：原始单元格和以“四角”点（2x2中间）为中心的相同尺寸。对于任何一对点，适当的距离是这两个中较小的一个。注意，这个方法也给出了一个框架，在这个框架中你可以很容易地计算动量变化。

第三种可能的方法是将尺寸增加一倍（如上面的2x2)，但给每个粒子四组坐标，每个方框中一组。在计算距离时，修改你的算法以考虑所有四个因素。如果您有良好的矢量化包和并行性，这可能是首选的解决方案。