当前位置 博文首页 > DECHIN：分子动力学模拟之周期性边界处理
周期性边界是分子动力学模拟中常用的一种技术手段，不仅可以完整的概述完整的分子体系的特性，在一部分场景下还可以提升计算的效率，从作用上来看更像是一类的近似模型（假设有一个原子逃出这个周期性边界封装的盒子，一定会有另一个相同原子从相对的边界走进这个盒子）。
首先我们用简单的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的动态图，效果如下：
在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的使用，在循环操作下也是有相当的编译优化空间的。
本文从分子动力学模拟中的周期性边界处理角度出发，介绍了无符号整数和带符号整数的一些应用的技巧，使用这些格式转换的技术有可能在程序的性能优化中带来一定的效果。同时为了更加直观的展示分子模拟的效果，我们用animation展示了一个简单的动态图绘制的案例，其中还包含了多个子图的绘图技术。
本文首发链接为：https://www.cnblogs.com/dechinphy/p/box.html
作者ID：DechinPhy
更多原著文章请参考：https://www.cnblogs.com/dechinphy/
打赏专用链接：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
腾讯云专栏同步：https://cloud.tencent.com/developer/column/91958
