分子動力學模擬之周期性邊界處理
- 2021 年 6 月 29 日
- 筆記
- Molecular Dynamics
技術背景
周期性邊界是分子動力學模擬中常用的一種技術手段，不僅可以完整的概述完整的分子體系的特性，在一部分場景下還可以提升計算的效率，從作用上來看更像是一類的近似模型（假設有一個原子逃出這個周期性邊界封裝的盒子，一定會有另一個相同原子從相對的邊界走進這個盒子）。
不加周期性邊界的場景
首先我們用簡單的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的使用，在循環操作下也是有相當的編譯優化空間的。
總結概要
本文從分子動力學模擬中的周期性邊界處理角度出發，介紹了無符號整數和帶符號整數的一些應用的技巧，使用這些格式轉換的技術有可能在程式的性能優化中帶來一定的效果。同時為了更加直觀的展示分子模擬的效果，我們用animation展示了一個簡單的動態圖繪製的案例，其中還包含了多個子圖的繪圖技術。
版權聲明
本文首發鏈接為：//www.cnblogs.com/dechinphy/p/box.html
作者ID：DechinPhy
更多原著文章請參考：//www.cnblogs.com/dechinphy/
打賞專用鏈接：//www.cnblogs.com/dechinphy/gallery/image/379634.html
騰訊雲專欄同步：//cloud.tencent.com/developer/column/91958
