分子动力学模拟之周期性边界处理

DechinPhy

发布于 2021-07-01 18:00:00

1.1K0

发布于 2021-07-01 18:00:00

技术背景

周期性边界是分子动力学模拟中常用的一种技术手段，不仅可以完整的概述完整的分子体系的特性，在一部分场景下还可以提升计算的效率，从作用上来看更像是一类的近似模型（假设有一个原子逃出这个周期性边界封装的盒子，一定会有另一个相同原子从相对的边界走进这个盒子）。

不加周期性边界的场景

首先我们用简单的python代码演示一个没加周期性边界条件的示例，一个红色的原子从坐标轴的0位置处移动到100的位置，但是盒子大小仅仅设置为20，这个大小也是我们的可见范围。也就是说，超过20之后我们就看不见这个原子了，具体代码实现如下所示：

import numpy as np
import time
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display

fig = plt.figure()
ims = []
for i in range(100):
    plt.xlim(0.0,20.0,5.0)
    x1 = i
    y1 = 5
    im = plt.plot(x1, y1, 'o', color="red")
    ims.append(im)
    time.sleep(0.1)

ani = animation.ArtistAnimation(fig, ims, repeat_delay=0)
ani.save('mol.gif', writer='pillow')

运行完成后会在当前路径下生成一个名为mol.gif的动态图，效果如下：

使用uint类型实现周期性边界

在python中可以用numpy的数据类型来转换给定的数据，而且性能有一定的保障。如果是用c++来编码，我们知道格式转换和移位操作之类的性能非常高，相比于数据的乘加运算而言，这种操作速度要快上许多。这里我们使用无符号的整型变量来处理周期性边界问题，我们用numpy的一些具体操作来看下无符号整数变量的一些对应操作：

In [1]: import numpy as np

In [2]: np.uint(2**16) # 默认的uint是32位
Out[2]: 65536

In [3]: np.uint(2**16+1) # 默认的uint是32位
Out[3]: 65537

In [4]: np.uint16(2**16) # 指定16位，超过最大数2**16-1之后归零
Out[4]: 0

In [5]: np.uint16(2**16-1) # 指定16位，不超过最大数2**16-1结果不变
Out[5]: 65535

In [6]: np.uint16(-1) # 下限为0，超过下限后从最大数2**16-1开始计算
Out[6]: 65535

In [7]: np.uint16(-2**16) # 越过一个0之后又达到了边界的0
Out[7]: 0

In [8]: np.int16(2**15-1) # 带符号整数的最大数是2**15-1，比无符号整数位少了一个比特位
Out[8]: 32767

In [9]: np.int16(2**15) # 超过最大数是从最小数-2**15开始计数
Out[9]: -32768

In [10]: np.int16(-2**15) # 不越过最小数结果不变
Out[10]: -32768

In [11]: np.int16(-2**15-1) # 越过最小数从最大数2**15-1开始计数
Out[11]: 32767

再回过头来思考一下其中的逻辑，首先，int16的一个比特位被用来做符号存储，因此最大可表示的数字是

2^{15}-1

，最小可表示的数字为

-2^{15}

。关于为什么负数的数量比正数多一个，这是因为16个比特位一共可以表示

2^{16}

个数字，那么如果包括0在内的话，只有在区分正0和负0的情况下，正数和负数的数量才会是一样的。所以，我们可以将负0可以用来表示

-2^{15}

这个数，这样看起来就多出来了一个负数，实际上只是一种优化的策略。无符号整数和带符号的整数都是周期性的锯齿形函数，但是无符号整数取得的空间都在正数上，所以在分子动力学模拟中更倾向于取无符号整数来处理周期性边界问题。

为了更加清晰的展现无符号整数的函数图像与周期性边界条件下的原子运动轨迹，我们将两张图画在一起来看下这个结果：

import numpy as np
import time
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display

fig = plt.figure()
ims = []
box_size = 20.0
x = []
y = []
for i in range(100):
    plt.subplot(211)
    x.append(i)
    y.append(np.uint16(i*(2**16-1)/box_size)*box_size/(2**16-1))
    plt.plot(x, y, color="black")
    
    im = plt.plot(x[-1], y[-1], 'o', color="red")
    ims.append(im)

    plt.subplot(212)
    plt.xlim(0.0,box_size,5.0)
    x1 = np.uint16(i*(2**16-1)/box_size)*box_size/(2**16-1)
    y1 = 5
    im = plt.plot(x1, y1, 'o', color="red")
    ims.append(im)
    time.sleep(0.1)

ani = animation.ArtistAnimation(fig, ims, repeat_delay=0)
ani.save('mol.gif', writer='pillow')

运行后生成的图片如下图所示：

需要注意的是，这里做类型转换之前，要将周期性盒子的边长转化到跟无符号整数位长度一致，才能够使用无符号整数来处理周期性边界问题，所以先后有两次单位转换。但是如果我们只是需要判断是否超出了边界，那就不需要做第二次的单位转换。值得一提的是，如果采用格式转换的形式来做计算，而免去if的使用，在循环操作下也是有相当的编译优化空间的。